CN105954803A - Post-stack seismic inversion method and device - Google Patents

Abstract

The invention provides a post-stack seismic inversion method and a post-stack seismic inversion device, wherein the method comprises the following steps: obtaining synthetic seismic data according to the initial model of the wave impedance and the actual seismic data; optimizing an initial model of wave impedance according to the actual seismic data and the synthetic seismic data to obtain first optimal wave impedance; carrying out high-pass filtering on the first optimal wave impedance to obtain a high-frequency inversion result of the wave impedance; obtaining an actual seismic data envelope according to the actual seismic data, and obtaining a synthetic seismic data envelope according to the synthetic seismic data; according to the initial model of the actual seismic data envelope and the synthetic seismic data envelope optimization wave impedance, obtaining a second optimal wave impedance; low-pass filtering is carried out on the second optimal wave impedance to obtain a low-frequency inversion result of the wave impedance; and obtaining a final inversion result of the wave impedance according to the high-frequency inversion result of the wave impedance and the low-frequency inversion result of the wave impedance. By implementing the method, the inversion accuracy is improved, and quantitative description is realized in reservoir prediction.

Description

Post-stack seismic inversion method and device

technical field

The present invention relates to the technical field of geophysical exploration, in particular to a post-stack seismic inversion method and device.

Background technique

Seismic exploration is to excite an artificial seismic source on the surface, and the vibration caused by the seismic source propagates underground in the form of seismic waves, and under certain conditions, it is reflected upwards and sent back to the surface, and then the seismic waves are reflected back by the geophone data to obtain seismic records (also called earthquakes data). Seismic data is obtained through the underground medium, which is bound to be affected by the physical properties of the underground medium (such as lithology, density, porosity, etc.), which includes not only reflected waves, but also noise, direct waves, surface waves, etc. To use seismic data to predict oil and gas in underground reservoirs, data processing is first required to remove various noises, direct waves and surface waves. The surface wave energy in seismic data is very strong, and the frequency is only within a dozen Hz, and it is usually removed in the frequency domain. The frequency of the processed seismic data is generally from a dozen hertz to tens of hertz, and the frequency value below a dozen hertz is very small.

Seismic inversion is an important topic in the field of reservoir geophysical exploration. Post-stack seismic inversion is the most commonly used technology for reservoir prediction. Usually, seismic data is converted into wave impedance profiles for reservoir prediction. Due to the small value of the seismic data below a dozen Hz, usually only the wave impedance data matching the seismic frequency band can be obtained in the post-stack seismic inversion, and the background value (that is, the low frequency value) of the subsurface wave impedance parameter cannot be obtained. , so it is difficult to quantitatively predict the reservoir. In the actual inversion, the usual practice is to use the logging data in the seismic work area for well seismic calibration, and under the constraints of the seismic horizon, use the calibrated logging data for mathematical interpolation to establish a low-frequency model of the physical parameters of the underground medium (also called the initial model), and then use the initial model for inversion to obtain relatively accurate physical parameters of the underground medium. However, no matter how advanced the interpolation method is, it is difficult to fully reflect the real spatial changes of the subsurface structure, especially in the exploration stage when the number of logging wells is scarce. Therefore, it is difficult to achieve quantitative description in reservoir prediction by using logging data interpolation to obtain the initial model for inversion.

Contents of the invention

The invention provides a post-stack seismic inversion method and device to solve the problem that it is difficult to realize quantitative description in reservoir prediction in the prior art.

In order to achieve the above object, an embodiment of the present invention provides a post-stack seismic inversion method, including: obtaining synthetic seismic data according to the initial model of wave impedance and actual seismic data; optimizing the initial model of wave impedance according to the actual seismic data and synthetic seismic data , to obtain the first optimal wave impedance; perform high-pass filtering on the first optimal wave impedance to obtain the high-frequency inversion result of wave impedance; obtain the actual seismic data envelope according to the actual seismic data, and obtain the synthetic seismic data according to the synthetic seismic data Envelope; according to the actual seismic data envelope and the synthetic seismic data envelope, optimize the initial model of wave impedance to obtain the second optimal wave impedance; perform low-pass filtering on the second optimal wave impedance to obtain the low-frequency inversion result of wave impedance ; Obtain the final inversion result of wave impedance according to the high-frequency inversion result of wave impedance and the low-frequency inversion result of wave impedance.

In one of the embodiments, the synthetic seismic data is obtained according to the initial model of wave impedance and the actual seismic data, including: extracting the corresponding wavelet in the target interval of the actual seismic data; obtaining the reflection coefficient according to the initial model; waves to obtain synthetic seismic data.

In one of the embodiments, the synthetic seismic data is obtained according to the reflection coefficient and the wavelet according to the following formula: d=R*w; wherein, the vector d is the synthetic seismic data; the vector R is the reflection coefficient, wherein, the first of the reflection coefficient R The i sampling point data R _i (i=1,2,3,...,n-1) is as follows: The last sampling point data R _n =R _n-1 ; Z _i is the data of the ith sampling point of the initial model of wave impedance; Z _i+1 is the data of the i+1 sampling point of the initial model of wave impedance; the vector w is The wavelet corresponding to the target interval of the actual seismic data.

In one of the embodiments, the initial wave impedance model is optimized according to the actual seismic data and the synthetic seismic data to obtain the first optimal wave impedance, including: calculating the error between the synthetic seismic data and the actual seismic data; When the range is set, the initial model of wave impedance is the first optimal wave impedance; when the error is not within the preset range, the number of modifications of the initial model of wave impedance is set to j, j=1, and proceed as follows Iterative processing: use the nonlinear global optimization algorithm to modify the initial model of wave impedance for the jth time to obtain the jth modified model of wave impedance; obtain the jth modified synthetic earthquake based on the jth modified model of wave impedance and the actual seismic data data; calculate the error between the jth modified synthetic seismic data and the actual seismic data; when the error is within the preset range, the jth modified model of wave impedance is the first optimal wave impedance; when the error is not in the preset When within the range, replace j in iteration processing with j+1.

In one of the embodiments, the error between the synthetic seismic data and the actual seismic data is calculated as follows: Among them, J _d is the error between the synthetic seismic data and the actual seismic data;

The vector d ₀ is the actual seismic data, as follows; d ₀ =[d ₀₁ d ₀₂ d ₀₃ ... d _0n ] ^T ; where, d _0i (i=1,2,3...n) is the actual seismic data d _0th Data of i sampling points; vector d is synthetic seismic data, as follows: d=[d ₁ d ₂ d ₃ …d _n ] ^T ; where, d _i (i=1,2,3...n) is synthetic seismic data The data of the i-th sampling point of data d.

In one of the embodiments, the actual seismic data envelope is obtained according to the actual seismic data in the following manner: e ₀ =[e ₀₁ e ₀₂ e ₀₃ ... e _0n ] ^T , Among them, the vector e ₀ is the actual seismic data envelope, e _0i (i=1,2,3...n) is the i-th sampling point data of the actual seismic data envelope e ₀ ; the Hilbert transform of d _0i as follows: The synthetic seismic data envelope is obtained from the synthetic seismic data as follows: e=[e ₁ e ₂ e ₃ ... e _n ] ^T , Among them, the vector e is the synthetic seismic data envelope, e _i (i=1,2,3...n) is the ith sample point data of the synthetic seismic data envelope e; the Hilbert transform of d _i as follows:

In one of the embodiments, the initial wave impedance model is optimized according to the actual seismic data envelope and the synthetic seismic data envelope, and the second optimal wave impedance is obtained, including: calculating the difference between the synthetic seismic data envelope and the actual seismic data envelope The error between; when the error is within the preset range, the initial model of wave impedance is the second optimal wave impedance; when the error is not within the preset range, the number of revisions of the initial model of wave impedance is set to k, k = 1, and perform iterative processing as follows: use the nonlinear global optimization algorithm to modify the initial model of wave impedance for the kth time, and obtain the kth modified model of wave impedance; according to the kth modified model of wave impedance and the actual earthquake The kth modified synthetic seismic data envelope is obtained from the data; the error between the kth modified synthetic seismic data envelope and the actual seismic data envelope is calculated; when the error is within the preset range, the kth modified model of wave impedance is the second optimal wave impedance; when the error is not within the preset range, replace k in the iterative process with k+1.

In one of the embodiments, the error between the synthetic seismic data envelope and the actual seismic data envelope is calculated as follows: Among them, J _e is the error between the synthetic seismic data envelope and the actual seismic data envelope.

In one of the embodiments, the final inversion result of the wave impedance is obtained according to the high-frequency inversion result of the wave impedance and the low-frequency inversion result of the wave impedance in the following manner: Z=Z _1h +Z _2l ; wherein, the vector Z is The final inversion result of wave impedance, the vector Z _1h is the high-frequency inversion result of wave impedance, and the vector Z _2l is the low-frequency inversion result of wave impedance.

An embodiment of the present invention also provides a post-stack seismic inversion device, including: a seismic data synthesis module, used to obtain synthetic seismic data according to the initial wave impedance model and actual seismic data; a high-frequency inversion module, used to Data and synthetic seismic data optimize the initial model of wave impedance to obtain the first optimal wave impedance; perform high-pass filtering on the first optimal wave impedance to obtain the high-frequency inversion result of wave impedance; the low-frequency inversion module is used to The actual seismic data envelope is obtained from the seismic data, and the synthetic seismic data envelope is obtained according to the synthetic seismic data; the initial model of wave impedance is optimized according to the actual seismic data envelope and the synthetic seismic data envelope, and the second optimal wave impedance is obtained; 2. Perform low-pass filtering on the optimal wave impedance to obtain the low-frequency inversion result of wave impedance; the post-stack seismic inversion module is used to obtain the final wave impedance according to the high-frequency inversion result of wave impedance and the low-frequency inversion result of wave impedance Inversion results.

In one embodiment, the seismic data synthesis module is specifically used to: extract corresponding wavelets in the actual seismic data target interval; obtain reflection coefficients according to the initial model; obtain synthetic seismic data according to the reflection coefficients and wavelets.

In one of the embodiments, the high-frequency inversion module is specifically used to obtain the first optimal wave impedance in the following manner: calculate the error between the synthetic seismic data and the actual seismic data; when the error is within a preset range, the wave impedance The initial model of impedance is the first optimal wave impedance; when the error is not within the preset range, the number of modifications of the initial model of wave impedance is set to j, j=1, and iterative processing is carried out as follows: using nonlinear The global optimization algorithm modifies the initial model of wave impedance for the jth time to obtain the jth modified model of wave impedance; obtain the jth modified synthetic seismic data according to the jth modified model of wave impedance and the actual seismic data; calculate the jth modified Modify the error between the synthetic seismic data and the actual seismic data; when the error is within the preset range, the jth modified model of wave impedance is the first optimal wave impedance; when the error is not within the preset range, iteratively j in processing is replaced by j+1.

In one of the embodiments, the low-frequency inversion module is specifically used to obtain the second optimal wave impedance in the following manner: obtain the actual seismic data envelope according to the actual seismic data, and obtain the synthetic seismic data envelope according to the synthetic seismic data; calculate The error between the synthetic seismic data envelope and the actual seismic data envelope; when the error is within the preset range, the initial wave impedance model is the second optimal wave impedance; when the error is not within the preset range, set The number of revisions of the initial model of wave impedance is k, k=1, and iterative processing is carried out in the following manner: the initial model of wave impedance is modified for the kth time using a nonlinear global optimization algorithm to obtain the kth modified model of wave impedance; According to the kth modified model of wave impedance and the actual seismic data, the kth modified synthetic seismic data envelope is obtained; the error between the kth modified synthetic seismic data envelope and the actual seismic data envelope is calculated; when the error is in the preset When the error is within the range, the kth modification model of the wave impedance is the second optimal wave impedance; when the error is not within the preset range, replace k in the iterative process with k+1.

In one of the embodiments, the post-stack seismic inversion module is specifically configured to obtain the final inversion result of wave impedance according to the high-frequency inversion result of wave impedance and the low-frequency inversion result of wave impedance in the following manner: Z=Z _1h +Z _2l ; where, the vector Z is the final inversion result of wave impedance, the vector Z _1h is the high-frequency inversion result of wave impedance, and the vector Z _2l is the low-frequency inversion result of wave impedance.

In the technical solution provided by the embodiment of the present invention, the inversion is first performed based on the initial model, and the high-frequency inversion result matching the seismic frequency band is obtained by high-pass filtering the inversion result, and then the envelope of the seismic data is calculated, and the initial model is used again to perform Inversion, the inversion results are low-pass filtered to obtain low-frequency inversion results, and the high-frequency inversion results and low-frequency inversion results are combined to form the final inversion result, which highlights the low-frequency magnitude of the seismic data below a dozen Hz and improves the inversion Accuracy, to achieve quantitative description in reservoir prediction.

Description of drawings

In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following will briefly introduce the accompanying drawings that need to be used in the description of the embodiments. Obviously, the accompanying drawings in the following description are only of the present invention. For some embodiments, those of ordinary skill in the art can also obtain other drawings based on these drawings without any creative effort.

Fig. 1 is the flowchart of post-stack seismic inversion method in the embodiment of the present invention;

Fig. 2 is the specific flowchart of step 101 in the embodiment of the present invention;

Fig. 3 is the specific flowchart of step 102 in the embodiment of the present invention;

Fig. 4 is the specific flowchart of step 105 in the embodiment of the present invention;

Fig. 5 is a schematic diagram of the actual wave impedance data of the inspection well used in the embodiment of the present invention;

Fig. 6 is a schematic diagram of comparison between the low-frequency part of the actual wave impedance data of the inspection well used in the embodiment of the present invention and the initial model of wave impedance;

Fig. 7 is a schematic diagram of the high-frequency part of the actual wave impedance data of the inspection well used in the embodiment of the present invention;

Fig. 8 is a schematic diagram of an actual seismic data corresponding to the coordinates of the inspection well used in the embodiment of the present invention;

Fig. 9 is the frequency spectrum of a piece of actual seismic data corresponding to the coordinates of the inspection well used in the embodiment of the present invention;

Fig. 10 is a schematic diagram of wavelets corresponding to actual seismic data in an embodiment of the present invention;

Fig. 11 is a schematic diagram of comparison between the first optimal wave impedance data and the actual wave impedance data of the inspection well used in the embodiment of the present invention;

Fig. 12 is a schematic diagram of the comparison between the high-frequency part of the first optimal wave impedance data and the high-frequency part of the actual wave impedance data in the embodiment of the present invention;

Fig. 13 is a schematic diagram of comparing the envelope of the actual seismic data calculated in the embodiment of the present invention with a piece of actual seismic data corresponding to the coordinates of the inspection well used;

Fig. 14 is a comparison schematic diagram of the frequency spectrum of the envelope of the actual seismic data calculated in the embodiment of the present invention and the frequency spectrum of a piece of actual seismic data corresponding to the coordinates of the used inspection well;

Fig. 15 is a schematic diagram of the comparison between the second optimal wave impedance data and the actual wave impedance data of the inspection well used in the embodiment of the present invention;

Fig. 16 is a schematic diagram of the comparison between the low-frequency part of the second optimal wave impedance data and the low-frequency part of the actual wave impedance data of the inspection well used in the embodiment of the present invention;

Fig. 17 is a schematic diagram of the comparison between the final inversion result of wave impedance and the actual wave impedance data of the inspection well used in the embodiment of the present invention;

Fig. 18 is a schematic structural diagram of a post-stack seismic inversion device in an embodiment of the present invention.

detailed description

The following will clearly and completely describe the technical solutions in the embodiments of the present invention with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only some, not all, embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without making creative efforts belong to the protection scope of the present invention.

The invention provides a post-stack seismic inversion method to solve the problem in the prior art that it is difficult to achieve quantitative description in reservoir prediction by using logging data interpolation to obtain an initial model for inversion.

Fig. 1 is a flowchart of a post-stack seismic inversion method in an embodiment of the present invention. As shown in Figure 1, post-stack seismic inversion methods include:

Step 101: Obtain synthetic seismic data according to the initial model of wave impedance and actual seismic data.

Step 102: Optimizing the initial model of wave impedance according to the actual seismic data and the synthetic seismic data to obtain the first optimal wave impedance.

Step 103: Perform high-pass filtering on the first optimal wave impedance to obtain a high-frequency inversion result of wave impedance.

Step 104: Obtain an envelope of actual seismic data according to the actual seismic data, and obtain an envelope of synthetic seismic data according to the synthetic seismic data.

Step 105: Optimizing the initial wave impedance model according to the actual seismic data envelope and the synthetic seismic data envelope to obtain a second optimal wave impedance.

Step 106: Perform low-pass filtering on the second optimal wave impedance to obtain a low-frequency inversion result of the wave impedance.

Step 107: Obtain the final inversion result of the wave impedance according to the high-frequency inversion result of the wave impedance and the low-frequency inversion result of the wave impedance.

In the specific implementation, the initial model of the wave impedance is established first. There are many ways to establish the initial model of wave impedance. For example, it is possible to collect the actual seismic data processed with amplitude preservation and the corresponding logging data in the seismic work area, carry out well seismic calibration according to the logging data, and carry out calibration on the calibrated logging data. interpolation.

A post-stack seismic inversion method provided by an embodiment of the present invention will be described in detail below in combination with specific embodiments, and an inspection well and its corresponding seismic data in the Xujiahe Formation in Tongnan, Sichuan Province will be taken as an example for illustration.

The initial model of the wave impedance is established as follows:

Acquire amplitude-preserved seismic data and corresponding logging data in the seismic work area, perform well-seismic calibration based on the known prior geological horizons and drilling horizons, and calibrate the acoustic time difference curve of the depth domain logging to the time domain, Interpolation is performed using calibrated log data to create an initial model of wave impedance.

The initial model of wave impedance can be represented by the vector Z ₀ :

Z ₀ =[Z ₁ Z ₂ Z ₃ ... Z _n ] ^T ,

Wherein, Z _i (i=1, 2, 3, . . . , n) is the i-th sampling point data, and the vector Z ₀ has n elements in total.

Fig. 5 is a schematic diagram of actual wave impedance data of inspection wells used in the embodiment of the present invention. As shown in Figure 5, the solid line in Figure 5 is the actual wave impedance data of the inspection well used in the present invention, and the data has a total of 192 sampling points. Fig. 6 is a schematic diagram of the comparison between the low frequency part of the actual acoustic impedance data of the inspection well used in the embodiment of the present invention and the initial model of acoustic impedance. As shown in Fig. 6, the solid line in Fig. 6 is the low-frequency part of the actual wave impedance data in Fig. 5 after low-pass filtering of 0 Hz-10 Hz, and the dotted line is the 0 Hz-10 Hz established by interpolation using other logging data in the seismic work area The wave impedance initial model Z ₀ . Fig. 7 is a schematic diagram of the high-frequency part of the actual wave impedance data of the inspection well used in the embodiment of the present invention. As shown in Fig. 7, the solid line in Fig. 7 is the high frequency part obtained after subtracting the low frequency part (solid line) in Fig. 6 from the actual wave impedance data in Fig. 5 . Comparing Figure 6, it can be seen that there is a large gap between the initial wave impedance model established by interpolation and the low frequency part of the actual wave impedance data.

There are many ways to obtain synthetic seismic data. For example, synthetic seismic data can be obtained as follows:

Fig. 2 is the specific flowchart of step 101 in the embodiment of the present invention; As shown in Fig. 2, step 101 specifically comprises:

Step 201: Extract corresponding wavelets in the target layer of actual seismic data.

Step 202: Obtain the reflection coefficient according to the initial model.

Step 203: Obtain synthetic seismic data according to the reflection coefficient and wavelet.

During specific implementation, first execute step 201; actual seismic data can be represented by vector d ₀ , as follows:

d ₀ =[d ₀₁ d ₀₂ d ₀₃ ... d _0n ] ^T ,

Wherein, d _0i (i=1, 2, 3...n) is the i-th sampling point data of the actual seismic data d ₀ , and the vector d ₀ has n elements in total.

The corresponding wavelet can be represented by a vector w.

Fig. 8 is a schematic diagram of an actual seismic data corresponding to the coordinates of the test well used in the embodiment of the present invention. As shown in Figure 8, the solid line in Fig. 8 is a piece of actual seismic data at the corresponding place of the inspection well coordinates used in the present invention, and the data has 192 sampling points in total; Fig. 9 is the corresponding place of the inspection well coordinates used in the embodiment of the present invention Frequency spectrum of an actual seismic data. As shown in Figure 9, the frequency spectrum is missing frequencies below 10 Hz. Fig. 10 is a schematic diagram of wavelets corresponding to actual seismic data in an embodiment of the present invention.

Then, step 202 is executed to obtain the reflection coefficient in various ways, for example, the reflection coefficient can be obtained as follows:

The reflection coefficient is calculated from the initial model, and the reflection coefficient can be represented by a vector R, where the i-th sampling point data R _i (i=1,2,3,...,n-1) of the reflection coefficient R is expressed as follows :

${\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; Z</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; Z</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}{{\mathrm{<span\; class="notranslate">\; Z</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; Z</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}}<span\; class="notranslate">\; ,</span>$

The last sampling point data R _n =R _n-1 ;

Among them, Z _i is the i-th sampling point data of the initial wave impedance model; Z _i+1 is the i+1-th sampling point data of the initial wave impedance model; R _n is the n-th sampling point data of the reflection coefficient.

Finally, execute step 203 to obtain synthetic data; during specific implementation, the synthetic seismic data can be obtained according to the following formula:

d=R*w;

The symbol "*" in the above formula is the convolution symbol, and the synthetic seismic data is the convolution of reflection coefficient and wavelet.

Among them, the vector d is the synthetic seismic data, as follows:

d=[d ₁ d ₂ d ₃ ... d _n ] ^T ;

d _i (i=1,2,3...n) is the i-th sampling point data of the synthetic seismic data d, and the vector d has n elements in total.

After obtaining the synthetic seismic data, step 102 is performed; Fig. 3 is a specific flowchart of step 102 in an embodiment of the present invention; as shown in Fig. 3 , step 102 specifically includes:

Step 301: Calculate the error between the synthetic seismic data and the actual seismic data.

Step 302: Determine whether the error is within a preset range.

Step 303: When the error is within the preset range, the initial wave impedance model is the first optimal wave impedance.

Step 304: When the error is not within the preset range, set the number of times of modification of the initial wave impedance model to j.

Step 305: Set j=1.

Step 306: Using the nonlinear global optimization algorithm to modify the initial model of the wave impedance for the jth time to obtain the jth modified model of the wave impedance.

Step 307: Obtain the jth modified synthetic seismic data according to the jth modified model of wave impedance and the actual seismic data.

Step 308: Calculate the error between the jth modified synthetic seismic data and the actual seismic data.

Step 309: Determine whether the error is within a preset range.

Step 310: When the error is within the preset range, the j-th modified model of wave impedance is the first optimal wave impedance.

Step 311 : when the error is not within the preset range, j=j+1, return to step 306 .

During specific implementation, step 301 is executed first, and the error between the synthetic seismic data and the actual seismic data can be obtained in the following manner:

${\mathrm{<span\; class="notranslate">\; J</span>}}_{\mathrm{<span\; class="notranslate">\; d</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; d</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{<span\; class="notranslate">\; \&Sigma;</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span>$

Among them, J _d is the error between the synthetic seismic data and the actual seismic data.

In an embodiment, the first optimal wave impedance can be represented by a vector Z ₁ , which has the same number of elements as the initial wave impedance model Z ₀ . Fig. 11 is a schematic diagram of the comparison between the first optimal wave impedance data and the actual wave impedance data of the test well used in the embodiment of the present invention. As shown in Figure 11, the dotted line in Figure 11 is the data of the _first optimal wave impedance Z obtained through inversion in this embodiment, and the solid line is the actual wave impedance of the inspection well used in the embodiment of the present invention in Figure 5 Data; comparison shows that the first optimal wave impedance and the actual wave impedance data still have large errors in some parts.

After obtaining the first optimal wave impedance, step 103 is performed; in an embodiment, the high-frequency inversion result (high-frequency part) of the wave impedance can be represented by a vector Z _1h , which has the same value as the first optimal wave impedance Z ₁ the same number of elements. Fig. 12 is a schematic diagram of the comparison between the high-frequency part of the first optimal wave impedance data and the high-frequency part of the actual wave impedance data in the embodiment of the present invention. As shown in Figure 12, the dotted line in Figure 12 is the high-frequency inversion result Z _1h obtained after the first optimal wave impedance Z ₁ is subjected to high-pass filtering, and the solid line is the actual wave impedance of the inspection well used in the present invention in Figure 7 The high-frequency part of the data; the comparison shows that the error between the high-frequency part of the first optimal wave impedance and the high-frequency part of the actual wave impedance data is very small.

After obtaining the high-frequency inversion result of wave impedance, execute step 104; during specific implementation, the actual seismic data envelope is obtained according to the actual seismic data in the following manner:

e ₀ =[e ₀₁ e ₀₂ e ₀₃ ... e _0n ] ^T ,

${\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>\sqrt{{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}\mathrm{<span\; class="notranslate">\; i</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; d</span>}}}_{\mathrm{<span\; class="notranslate">\; 0</span>}\mathrm{<span\; class="notranslate">\; i</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; ;</span>$

Wherein, vector e ₀ is the actual seismic data envelope, and e _0i (i=1,2,3...n) is the i-th sampling point data of the actual seismic data envelope e ₀ ;

Hilbert transform of d _0i as follows:

${\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; d</span>}}}_{\mathrm{<span\; class="notranslate">\; 0</span>}\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; *</span>{\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; ,</span>$

${\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; -</span>{<span\; class="notranslate">\; (</span><span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; i</span>}}}{\mathrm{<span\; class="notranslate">\; i</span>}\mathrm{<span\; class="notranslate">\; \&pi;</span>}}<span\; class="notranslate">\; .</span>$

During specific implementation, the synthetic seismic data envelope is obtained according to the synthetic seismic data in the following manner:

e=[e ₁ e ₂ e ₃ ... e _n ] ^T ,

${\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>\sqrt{{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>{\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; d</span>}}}_{\mathrm{<span\; class="notranslate">\; i</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}<span\; class="notranslate">\; ;</span>$

Wherein, the vector e is the synthetic seismic data envelope, e _i (i=1,2,3...n) is the ith sampling point data of the synthetic seismic data envelope e;

Hilbert transform of d _i as follows:

${\stackrel{<span\; class="notranslate">\; ^</span>}{\mathrm{<span\; class="notranslate">\; d</span>}}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; d</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; *</span>{\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; .</span>$

Fig. 13 is a schematic diagram of the comparison between the envelope of the actual seismic data calculated in the embodiment of the present invention and a piece of actual seismic data corresponding to the coordinates of the inspection well used. As shown in Figure 13, the dotted line in Figure 13 is the envelope e ₀ of the actual seismic data calculated in the present invention, and the solid line is the actual seismic data corresponding to the coordinates of the inspection well in Figure 8; it can be seen from the comparison that the envelope e 0 The network can well reflect the reflection characteristics and trends of actual seismic data. Fig. 14 is a schematic diagram of the comparison between the frequency spectrum of the envelope of the actual seismic data calculated in the embodiment of the present invention and the frequency spectrum of a piece of actual seismic data corresponding to the coordinates of the inspection well used. As shown in Figure 14, the dotted line in Figure 14 is the frequency spectrum of the envelope _e0 of the actual seismic data calculated in the present invention, and the solid line is the frequency spectrum of the actual seismic data in Figure 9; as can be seen from the comparison, the envelope The frequency spectrum of the network well compensates for the low frequency part of the frequency spectrum of the actual seismic data.

After obtaining the actual seismic data envelope and the synthetic seismic data envelope, step 105 is performed; Fig. 4 is a schematic diagram of the actual wave impedance data of the test well used in the embodiment of the present invention; as shown in Fig. 4, step 105 specifically includes:

Step 401: Calculate the error between the synthetic seismic data envelope and the actual seismic data envelope.

Step 402: Determine whether the error is within a preset range.

Step 403: When the error is within the preset range, the initial wave impedance model is the second optimal wave impedance.

Step 404: When the error is not within the preset range, set the number of times of modification of the initial wave impedance model to k.

Step 405: Set k=1.

Step 406: Using the nonlinear global optimization algorithm to modify the initial model of the wave impedance for the kth time to obtain the kth time modified model of the wave impedance.

Step 407: According to the k-th modified model of wave impedance and the actual seismic data, the k-th modified synthetic seismic data envelope is obtained.

Step 408: Calculating the error between the kth modified synthetic seismic data envelope and the actual seismic data envelope.

Step 409: Determine whether the error is within a preset range.

Step 410: When the error is within the preset range, the kth modified model of wave impedance is the second optimal wave impedance.

Step 411 : when the error is not within the preset range, k=k+1, return to step 406 .

During specific implementation, step 401 is executed first, and the error between the synthetic seismic data envelope and the actual seismic data envelope can be obtained as follows:

${\mathrm{<span\; class="notranslate">\; J</span>}}_{\mathrm{<span\; class="notranslate">\; e</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&Sigma;</span>}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\underset{\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; no</span>}}{<span\; class="notranslate">\; \&Sigma;</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; e</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span>$

Among them, J _e is the error between the synthetic seismic data envelope and the actual seismic data envelope.

In an embodiment, the second optimal wave impedance can be represented by a vector Z ₂ , which has the same number of elements as the initial wave impedance model Z ₀ . Fig. 15 is a schematic diagram of the comparison between the second optimal wave impedance data and the actual wave impedance data of the test well used in the embodiment of the present invention. As shown in Figure 15, dotted line in Fig. 15 is the data of the _second optimum wave impedance Z2 that obtains through envelope inversion in the present invention, and solid line is the actual wave impedance data of inspection well used in Fig. 5; Comparison can be It can be seen that there are still large errors in some parts between the second optimal wave impedance and the actual wave impedance data.

After obtaining the second optimal wave impedance, step 106 is performed; in an embodiment, the low-frequency inversion result (low frequency part) of the wave impedance can be represented by a vector Z ₂₁ , which has the same value as the second optimal wave impedance Z ₂ number of elements. Fig. 16 is a schematic diagram of the comparison between the low-frequency part of the second optimal wave impedance data and the low-frequency part of the actual wave impedance data of the test well used in the embodiment of the present invention. As shown in Fig. 16, the dotted line in Fig. 16 is the low-frequency inversion result Z _2l obtained after the second optimal wave impedance Z ₂ is low-pass filtered, and the solid line is the low-frequency of the actual wave impedance data of the test well used in Fig. 6 part (solid line); it can be seen from the comparison that the error between the low frequency part of the second optimal wave impedance and the low frequency part of the actual wave impedance data is very small.

After obtaining the high-frequency inversion result of wave impedance, execute step 107; wherein, the final inversion result of wave impedance can be represented by vector Z:

Z = Z _1h + Z _2l .

Fig. 17 is a schematic diagram of the comparison between the final inversion result of wave impedance and the actual wave impedance data of the inspection well used in the embodiment of the present invention. As shown in Fig. 17, the dotted line in Fig. 17 is the final inversion result Z of wave impedance, and the solid line is the actual wave impedance data of the inspection well used in Fig. The final inversion result of the low-frequency part of network inversion largely restores the real wave impedance.

Based on the same inventive concept, an embodiment of the present invention also provides a post-stack seismic inversion device. Since the problem-solving principle of the device is similar to that of the post-stack seismic inversion method, the implementation of the device can refer to the implementation of the method, repeat The place will not be repeated.

Fig. 18 is a schematic structural diagram of a post-stack seismic inversion device in an embodiment of the present invention. As shown in Fig. 18, the device may include:

Seismic data synthesis module 181, used for obtaining synthetic seismic data according to the initial model of wave impedance and actual seismic data;

The high-frequency inversion module 182 is used to optimize the initial model of wave impedance according to the actual seismic data and synthetic seismic data to obtain the first optimal wave impedance; perform high-pass filtering on the first optimal wave impedance to obtain the high-frequency inversion of wave impedance performance results;

The low-frequency inversion module 183 is used to obtain the actual seismic data envelope according to the actual seismic data, and obtain the synthetic seismic data envelope according to the synthetic seismic data; optimize the initial model of wave impedance according to the actual seismic data envelope and the synthetic seismic data envelope, obtaining the second optimal wave impedance; performing low-pass filtering on the second optimal wave impedance to obtain a low-frequency inversion result of the wave impedance;

The post-stack seismic inversion module 184 is configured to obtain the final inversion result of wave impedance according to the high-frequency inversion result of wave impedance and the low-frequency inversion result of wave impedance.

In a specific embodiment of the present invention, the seismic data synthesis module is specifically used to: extract the corresponding wavelet in the actual seismic data target layer; obtain the reflection coefficient according to the initial model; obtain the synthesis according to the reflection coefficient and the wavelet seismic data.

In a specific embodiment of the present invention, the high-frequency inversion module is specifically used to obtain the first optimal wave impedance in the following manner: calculate the error between the synthetic seismic data and the actual seismic data; when the error is within a preset range , the initial model of wave impedance is the first optimal wave impedance; when the error is not within the preset range, the number of revisions of the initial model of wave impedance is set to j, j=1, and iterative processing is carried out as follows: The nonlinear global optimization algorithm modifies the initial model of wave impedance for the jth time to obtain the jth modified model of wave impedance; obtain the jth modified synthetic seismic data according to the jth modified model of wave impedance and the actual seismic data; calculate the jth modified The error between the synthetic seismic data and the actual seismic data is modified j times; when the error is within the preset range, the jth modified model of wave impedance is the first optimal wave impedance; when the error is not within the preset range, Replace j in iterative processing with j+1.

In a specific embodiment of the present invention, the low-frequency inversion module is specifically used to obtain the second optimal wave impedance as follows:

The actual seismic data envelope is obtained from the actual seismic data, and the synthetic seismic data envelope is obtained from the synthetic seismic data;

calculating the error between the synthetic seismic data envelope and the actual seismic data envelope;

When the error is within the preset range, the initial wave impedance model is the second optimal wave impedance;

When the error is not within the preset range, set the number of revisions of the initial model of wave impedance to k, k=1, and perform iterative processing as follows:

Using the nonlinear global optimization algorithm to modify the initial model of wave impedance for the kth time, and obtain the kth modified model of wave impedance;

According to the kth modified model of wave impedance and the actual seismic data, the kth modified synthetic seismic data envelope is obtained;

Calculating the error between the kth modified synthetic seismic data envelope and the actual seismic data envelope;

When the error is within the preset range, the kth modified model of wave impedance is the second optimal wave impedance;

When the error is not within the preset range, k in the iterative process is replaced by k+1.

In a specific embodiment of the present invention, the post-stack seismic inversion module is specifically used to obtain the final inversion result of wave impedance according to the high-frequency inversion result of wave impedance and the low-frequency inversion result of wave impedance in the following manner: Z= Z _1h +Z _2l ; where, the vector Z is the final inversion result of wave impedance, the vector Z _1h is the high-frequency inversion result of wave impedance, and the vector Z _2l is the low-frequency inversion result of wave impedance.

It can be seen from the above-mentioned embodiments that in the technical solution provided by the embodiment of the present invention, the inversion is first performed based on the initial model, the inversion result is high-pass filtered to obtain the high-frequency inversion result matching the seismic frequency band, and then the envelope of the seismic data is calculated , using the initial model again for inversion, low-pass filtering the inversion results to obtain the low-frequency inversion results, combining the high-frequency inversion results and low-frequency inversion results into the final inversion result, highlighting the low-frequency content of the seismic data below a dozen Hz It solves the problem that the inversion results are difficult to quantify due to the low frequency and magnitude of seismic data, greatly improves the inversion accuracy, and has obvious effects in the prediction of complex underground reservoirs.

Those skilled in the art should understand that the embodiments of the present invention may be provided as methods, systems, or computer program products. Accordingly, the present invention can take the form of an entirely hardware embodiment, an entirely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, the present invention may take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) having computer-usable program code embodied therein.

The present invention is described with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It should be understood that each procedure and/or block in the flowchart and/or block diagram, and a combination of procedures and/or blocks in the flowchart and/or block diagram can be realized by computer program instructions. These computer program instructions may be provided to a general purpose computer, special purpose computer, embedded processor, or processor of other programmable data processing equipment to produce a machine such that the instructions executed by the processor of the computer or other programmable data processing equipment produce a An apparatus for realizing the functions specified in one or more procedures of the flowchart and/or one or more blocks of the block diagram.

These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing apparatus to operate in a specific manner, such that the instructions stored in the computer-readable memory produce an article of manufacture comprising instruction means, the instructions The device realizes the function specified in one or more procedures of the flowchart and/or one or more blocks of the block diagram.

These computer program instructions can also be loaded onto a computer or other programmable data processing device, causing a series of operational steps to be performed on the computer or other programmable device to produce a computer-implemented process, thereby The instructions provide steps for implementing the functions specified in the flow chart or blocks of the flowchart and/or the block or blocks of the block diagrams.

The specific embodiments described above have further described the purpose, technical solutions and beneficial effects of the present invention in detail. It should be understood that the above descriptions are only specific embodiments of the present invention and are not intended to limit the scope of the present invention. Protection scope, within the spirit and principles of the present invention, any modification, equivalent replacement, improvement, etc., shall be included in the protection scope of the present invention.

Claims (14)

1. a poststack seismic inversion method, it is characterised in that including:

Initial model according to natural impedance and actual seismic data obtain synthetic seismic data；

Optimize the initial model of natural impedance according to actual seismic data and synthetic seismic data, obtain the first optimum natural impedance；

First optimum natural impedance is carried out high-pass filtering, obtains the high frequency inversion result of natural impedance；

Obtain actual seismic data envelopment according to actual seismic data, and obtain synthetic seismic data bag according to synthetic seismic data Network；

According to the initial model of actual seismic data envelopment Yu synthetic seismic data envelop optimization natural impedance, obtain the second optimum ripple Impedance；

Second optimum natural impedance is carried out low-pass filtering, obtains the low frequency inversion result of natural impedance；

High frequency inversion result according to natural impedance and the low frequency inversion result of natural impedance obtain the final inversion result of natural impedance.

Poststack seismic inversion method the most according to claim 1, it is characterised in that the described initial model according to natural impedance Synthetic seismic data is obtained with actual seismic data, including:

Corresponding wavelet is extracted at actual seismic data interval of interest；

Reflection coefficient is obtained according to initial model；

Synthetic seismic data is obtained according to described reflection coefficient and described wavelet.

Poststack seismic inversion method the most according to claim 2, it is characterised in that as follows according to described reflection system Number and described wavelet obtain synthetic seismic data:

D=R*w；

Wherein, vector d is synthetic seismic data；

Vector R is reflection coefficient, wherein, and ith sample point data R of reflection R_i(i=1,2,3 ..., n-1) as follows:

$R_i=\frac{Z_{i+1}-Z_i}{Z_{i+1}+Z_i};$

Last sampling number is according to R_n=R_n-1；

Z_iInitial model ith sample point data for natural impedance；Z_i+1Initial model i+1 sampling number for natural impedance According to；

Vector w is the wavelet that actual seismic data interval of interest is corresponding.

Poststack seismic inversion method the most according to claim 1, it is characterised in that described according to actual seismic data and conjunction Become geological data to optimize the initial model of natural impedance, obtain the first optimum natural impedance, including:

Calculate the error between synthetic seismic data and actual seismic data；

When error is in preset range, the initial model of natural impedance is the first optimum natural impedance；

When error is not in preset range, set the amendment number of times of initial model of natural impedance as j, j=1, and according to such as Under type is iterated processing:

Use the initial model of non-linear global optimizing algorithm jth time amendment natural impedance, obtain the jth time amendment mould of natural impedance Type；

Jth according to natural impedance time amendment model and actual seismic data obtain jth time amendment synthetic seismic data；

Calculate the error between jth time amendment synthetic seismic data and actual seismic data；

When error is in preset range, the jth of natural impedance time amendment model is the first optimum natural impedance；

When error is not in preset range, the j in iterative processing is replaced by j+1.

Poststack seismic inversion method the most according to claim 4, it is characterised in that calculate synthesis earthquake as follows Error between data and actual seismic data:

$J_d=\mathrm{\&Sigma;}{(d_0-d)}^2=\underset{i=1}{\overset{n}{\&Sigma;}}{(d_{0i}-d_i)}^2;$

Wherein, J_dFor the error between synthetic seismic data and actual seismic data；

Vector d₀For actual seismic data, as follows；

d₀=[d₀₁ d₀₂ d₀₃ … d_0n]^T；

Wherein, d_0i(i=1,2,3...n) it is actual seismic data d₀Ith sample point data；

Vector d is synthetic seismic data, as follows:

D=[d₁ d₂ d₃ … d_n]^T；

Wherein, d_i(i=1,2,3...n) it is synthetic seismic data d ith sample point data.

Poststack seismic inversion method the most according to claim 5, it is characterised in that as follows according to actual seismic Data obtain actual seismic data envelopment:

e₀=[e₀₁ e₀₂ e₀₃ … e_0n]^T,

$e_{0i}=\sqrt{d_{0i}^2+{\hat{d}}_{0i}^2};$

Wherein, vector e₀For actual seismic data envelopment, e_0i(i=1,2,3...n) it is actual seismic data envelopment e₀I-th is adopted Sampling point data；

d_0iHilbert transformAs follows:

${\hat{d}}_{0i}=d_{0i}*h_i,$

$h_i=\frac{1-{(-1)}^i}{i\mathrm{\&pi;}};$

Synthetic seismic data envelope is obtained as follows according to synthetic seismic data:

E=[e₁ e₂ e₃ … e_n]^T,

$e_i=\sqrt{d_i^2+{\hat{d}}_i^2};$

Wherein, vector e is synthetic seismic data envelope, e_i(i=1,2,3...n) it is synthetic seismic data envelope e ith sample Point data；

d_iHilbert transformAs follows:

${\hat{d}}_i=d_i*h_i.$

Poststack seismic inversion method the most according to claim 6, it is characterised in that described according to actual seismic data envelopment With the initial model of synthetic seismic data envelop optimization natural impedance, obtain the second optimum natural impedance, including:

Calculate the error between synthetic seismic data envelope and actual seismic data envelopment；

When error is in preset range, the initial model of natural impedance is the second optimum natural impedance；

When error is not in preset range, set the amendment number of times of initial model of natural impedance as k, k=1, and according to such as Under type is iterated processing:

Use the initial model of non-linear global optimizing algorithm kth time amendment natural impedance, obtain the kth time amendment mould of natural impedance Type；

Kth according to natural impedance time amendment model and actual seismic data obtain kth time amendment synthetic seismic data envelope；

Calculate the error between kth time amendment synthetic seismic data envelope and actual seismic data envelopment；

When error is in preset range, the kth of natural impedance time amendment model is the second optimum natural impedance；

When error is not in preset range, the k in iterative processing is replaced by k+1.

Poststack seismic inversion method the most according to claim 7, it is characterised in that calculate synthesis earthquake as follows Error between data envelopment and actual seismic data envelopment:

$J_e=\mathrm{\&Sigma;}{(e_0-e)}^2=\underset{i=1}{\overset{n}{\&Sigma;}}{(e_{0i}-e_i)}^2;$

Wherein, J_eFor the error between synthetic seismic data envelope and actual seismic data envelopment.

Poststack seismic inversion method the most according to claim 1, it is characterised in that as follows according to natural impedance The low frequency inversion result of high frequency inversion result and natural impedance obtains the final inversion result of natural impedance:

Z=Z_1h+Z_2l；

Wherein, vector Z is the final inversion result of natural impedance, vector Z_1hFor the high frequency inversion result of natural impedance, vector Z_2lFor ripple The low frequency inversion result of impedance.

10. a poststack seismic inversion device, it is characterised in that including:

Geological data synthesis module, obtains synthetic seismic data for the initial model according to natural impedance and actual seismic data；

High frequency inverting module, for optimizing the initial model of natural impedance according to actual seismic data and synthetic seismic data, obtains First optimum natural impedance；First optimum natural impedance is carried out high-pass filtering, obtains the high frequency inversion result of natural impedance；

Low frequency inverting module, for obtaining actual seismic data envelopment according to actual seismic data, and according to synthetic seismic data Obtain synthetic seismic data envelope；Introductory die according to actual seismic data envelopment Yu synthetic seismic data envelop optimization natural impedance Type, obtains the second optimum natural impedance；Second optimum natural impedance is carried out low-pass filtering, obtains the low frequency inversion result of natural impedance；

Poststack seismic inversion module, for obtaining ripple according to the high frequency inversion result of natural impedance and the low frequency inversion result of natural impedance The final inversion result of impedance.

11. poststack seismic inversion devices according to claim 10, it is characterised in that described geological data synthesis module has Body is used for:

Corresponding wavelet is extracted at actual seismic data interval of interest；

Reflection coefficient is obtained according to initial model；

Synthetic seismic data is obtained according to described reflection coefficient and described wavelet.

12. poststack seismic inversion devices according to claim 10, it is characterised in that described high frequency inverting module is specifically used In obtaining the first optimum natural impedance as follows:

Calculate the error between synthetic seismic data and actual seismic data；

When error is in preset range, the initial model of natural impedance is the first optimum natural impedance；

When error is not in preset range, set the amendment number of times of initial model of natural impedance as j, j=1, and according to such as Under type is iterated processing:

Use the initial model of non-linear global optimizing algorithm jth time amendment natural impedance, obtain the jth time amendment mould of natural impedance Type；

Jth according to natural impedance time amendment model and actual seismic data obtain jth time amendment synthetic seismic data；

Calculate the error between jth time amendment synthetic seismic data and actual seismic data；

When error is in preset range, the jth of natural impedance time amendment model is the first optimum natural impedance；

When error is not in preset range, the j in iterative processing is replaced by j+1.

13. poststack seismic inversion devices according to claim 10, it is characterised in that described low frequency inverting module is specifically used In obtaining the second optimum natural impedance as follows:

Obtain actual seismic data envelopment according to actual seismic data, and obtain synthetic seismic data bag according to synthetic seismic data Network；

Calculate the error between synthetic seismic data envelope and actual seismic data envelopment；

When error is in preset range, the initial model of natural impedance is the second optimum natural impedance；

When error is not in preset range, set the amendment number of times of initial model of natural impedance as k, k=1, and according to such as Under type is iterated processing:

Use the initial model of non-linear global optimizing algorithm kth time amendment natural impedance, obtain the kth time amendment mould of natural impedance Type；

Kth according to natural impedance time amendment model and actual seismic data obtain kth time amendment synthetic seismic data envelope；

Calculate the error between kth time amendment synthetic seismic data envelope and actual seismic data envelopment；

When error is in preset range, the kth of natural impedance time amendment model is the second optimum natural impedance；

When error is not in preset range, the k in iterative processing is replaced by k+1.

14. poststack seismic inversion devices according to claim 10, it is characterised in that described poststack seismic inversion module has Body is for obtaining natural impedance according to the high frequency inversion result of natural impedance and the low frequency inversion result of natural impedance as follows Final inversion result:

Z=Z_1h+Z_2l；

Wherein, vector Z is the final inversion result of natural impedance, vector Z_1hFor the high frequency inversion result of natural impedance, vector Z_2lFor ripple The low frequency inversion result of impedance.