CN105388520A - Seismic data pre-stack reverse time migration imaging method - Google Patents

Abstract

The invention discloses a seismic data pre-stack reverse time migration imaging method, which belongs to the technical field of oil exploration and development. The seismic data pre-stack reverse time migration imaging method utilizes physical property parameters represented by a three-dimensional depth domain speed field, performs partitioning processing on a work area, determines grid scales of different areas, performs grid discrete processing, utilizes a GPU to share the massive computation of wave field extrapolation, utilizes the feature that reverse time migration and shot-domain migration are irrelevant and use a CPU to perform shot-domain coarse-grain parallelism, and effectively increases imaging precision and imaging efficiency of pre-stack reverse time migration.

Description

A kind of seismic data pre-Stack Reverse formation method

Technical field

The present invention discloses a kind of seismic data pre-Stack Reverse formation method, belongs to petroleum exploration and development technical field.

Background technology

Pre-Stack Reverse is based on all-wave field wave equation, to single shot record reverse extrapolation in time, realizes underground each point accurately image, improves the imaging precision of skew achievement.Pre-Stack Reverse is when utilizing the discrete grid block of rule to carry out wave field extrapolation, discrete grid block yardstick is less, pre-Stack Reverse imaging precision is higher, but the calculated amount that pre-Stack Reverse wave field calculates and memory amount will be larger, imaging efficiency is lower, calculated amount and memory amount even can exceed the amount of bearing of computing machine, cause pre-Stack Reverse work to carry out.Otherwise discrete grid block yardstick is larger, calculated amount and the memory amount of wave field calculating are less, and imaging efficiency is higher, but pre-Stack Reverse imaging precision will be lower, can not reflect underground looks well, cannot embody the advantage of pre-Stack Reverse.The contradiction of pre-Stack Reverse imaging precision and imaging efficiency has become the principal contradiction restricting pre-Stack Reverse and extensively use in oil field.

Along with the development of computer technology, graphic based processor GPU computing platform is accelerated reverse-time migration formation method, becomes the effective means promoting reverse-time migration imaging efficiency.But, the storage space of GPU own is less, and local storage maximumly only has 6GB, when calculated amount and memory space excessive time, still need to adopt the calculated amount and the memory amount that reduce pre-Stack Reverse compared with large scale grid, in order to ensure the smooth enforcement of pre-Stack Reverse.Therefore, when portraying crack, hole small scale geologic body, minor fault, imaging precision still has much room for improvement.In recent years, based on GPU computing platform, these problems be also reflects to some achievements that reverse-time migration formation method accelerates.

Patent " a kind of based on GPU little memory space staggered-mesh 3-D seismics pre-Stack Reverse formation method (CN201010188322) " discloses a kind of based on GPU little memory space staggered-mesh 3-D seismics pre-Stack Reverse formation method, utilize the parallel computing based on GPU hardware and the technology such as Difference Schemes with Staggered, checkpoint storage, bow-tie and wide-angle reflection ripple accurate imaging can be made.This technology mainly utilizes GPU to carry out little memory space enforcement, and the storage by the highest 6GB of GPU limits, and the model data less to data volume can only carry out pre-Stack Reverse, is difficult to the actual seismic data processing magnanimity.

Geophysical prospecting for oil magazine 48 volume the 5th phase discloses one " the reverse-time migration formation method based on GPU parallel accelerate ", the method establishes a set of pre-Stack Reverse method accelerated based on GPU, source reconstruction method is used to solve reverse-time migration algorithm mass data storage and the large problem of calculated amount, and by the method Successful utilization to Sichuan Basin block.The grid that GPU accelerates only, shape relatively large at yardstick rule relatively improves imaging efficiency, but portray crack, hole small scale geologic body and complicated minor fault time, its pre-Stack Reverse imaging precision is relatively low.

In sum, under existence conditions, for magnanimity 3D seismic data, particularly high accuracy surface modelling data, only rely on GPU hardware computing platform, still can not improve the imaging precision of pre-Stack Reverse, the contradiction between the imaging precision of pre-Stack Reverse and imaging efficiency still exists.

Summary of the invention

The present invention seeks to overcome above-mentioned pre-Stack Reverse imaging technique imaging precision and the relatively low defect of imaging efficiency, a kind of Moving grids pre-Stack Reverse method based on GPU parallel accelerate of proposition.

For achieving the above object, the present invention is by the following technical solutions:

1, utilize the three-dimensional earthquake data before superposition in work area to be evaluated, set up its three dimensional depth territory velocity field.

2, utilize three dimensional depth territory, work area to be evaluated velocity field, set up the Moving grids model of three dimensional depth territory velocity field different scale.

The physical parameter in the work area to be evaluated that the 2.1 three dimensional depth territory velocity fields utilizing step 1 to generate characterize, carries out subregion to work area to be evaluated, determines the mesh scale of different subregion;

2.2, by the discrete processes to the mesh scale of the different subregions that step 2.1 is determined, set up the Moving grids model of three dimensional depth territory velocity field.

3, set man-made explosion, just drilling its propagation property in the Moving grids model of three dimensional depth territory velocity field different scale, determining man-made explosion wave field, obtain its border wave field, and be saved in CPU.

The man-made explosion that 3.1 settings match with three-dimensional earthquake data before superposition dominant frequency, is just drilling its propagation property in the Moving grids model of three dimensional depth territory velocity field different scale, is determining man-made explosion wave field;

The border wave field of man-made explosion wave field in 3.2 obtaining steps 3.1, and be saved in CPU.

4, utilize three-dimensional earthquake data before superposition, the inverse time extrapolates its propagation property in the Moving grids model of three dimensional depth territory velocity field different scale, obtains the inverse time extrapolation wave field of three-dimensional earthquake data before superposition.

5, utilize step 3.2 to be kept at the border wave field of the man-made explosion wave field in CPU, on GPU, the inverse time of the three-dimensional earthquake data before superposition that reconstruct generates with step 4 extrapolates wave field man-made explosion wave field in the same time.

The inverse time extrapolation wave field of the three-dimensional earthquake data before superposition that the man-made explosion wave field 6, utilizing step 5 to reconstruct and step 4 generate, determines single big gun reverse-time migration imaging results, and then determines the reverse-time migration imaging results in work area to be evaluated.

The invention has the beneficial effects as follows: the present invention utilize three dimensional depth territory velocity field to characterize physical parameter, multidomain treat-ment is carried out to work area, determine the mesh scale of zones of different, and carry out Moving grids discrete processes, the magnanimity utilizing GPU to share wave field extrapolation calculates, utilizing reverse-time migration big gun territory to offset incoherent feature uses CPU to carry out big gun territory coarse grain parallelism, effectively improves imaging precision and the imaging efficiency of pre-Stack Reverse.

Accompanying drawing explanation

Fig. 1 is technical solution of the present invention FB(flow block);

Fig. 2 is the Moving grids model of three dimensional depth territory velocity field different scale;

Fig. 3 is Dong Pu oil field hollow shallow-layer small complex fault block district's reverse-time migration imaging results;

Fig. 4 is the reverse-time migration imaging results of Dong Pu oil field hollow mid-deep strata height steep fault;

Fig. 5 is hollow shallow-layer small complex fault block district, Dong Pu oil field depth migration imaging result;

Fig. 6 is the depth migration imaging result of Dong Pu oil field hollow mid-deep strata height steep fault.

Embodiment

With Zhongyuan Oil Field Branch, China Petrochemical Corp., Ltd. Dong Pu oil field hollow high accuracy surface modelling data imaging processing for example with reference, the present invention is described further.The embodiment of the present invention is as follows as shown in Figure 1:

1, obtain hollow work area, Dong Pu oil field high-precision three-dimensional earthquake data before superposition, its dominant frequency 30hz, frequency range is 3hz-65hz, and bin is 5m × 10m, and degree of covering is 8 × 4.Utilize three-dimensional earthquake data before superposition to set up its three dimensional depth territory velocity field, the maximal rate of three dimensional depth territory velocity field is 6100m/s, and minimum speed is 1900m/s, and the degree of depth is 13000m.

2, utilize three dimensional depth territory, work area velocity field, set up the Moving grids model of three dimensional depth territory velocity field different scale.

Work area speed, dominant frequency physical parameter that the 2.1 three dimensional depth territory velocity fields utilizing step 1 to generate characterize, carry out subregion to work area.According to the velocity variations feature of three dimensional depth territory velocity field along geologic horizon Dongying Formation, four groups, sand, certain fracture, work area is divided into the shallow-layer low regime on Dongying Formation, the small complex fault block district between Dongying Formation and substrate and three, the mid-deep strata district subregion under the steep major fault of hollow target area, substrate and height, each subregion minimum speed is: shallow-layer low regime minimum speed is 1900m/s, small complex fault block district and hollow target area minimum speed are 3000m/s, and mid-deep strata district minimum speed is 5600m/s.Formula (1) is utilized to determine that the mesh scale of different subregion is respectively: 6m, 10m, 18m.

$\mathrm{\&Delta;}x=\mathrm{\&Delta;}y<0.2*\frac{V_{min}}{2*feq}---\left(1\right)$

Wherein, Δ x is the discrete yardstick of grid in x direction, and Δ y is the discrete yardstick of grid in y direction, v _minfor different subregion minimum speed, feq is dominant frequency of seismic wavelet.

2.2, by the discrete processes to the mesh scale of three subregions that step 2.1 is determined, set up the Moving grids model of three dimensional depth territory velocity field three kinds of different scales from shallow to deep.As shown in Figure 2: I district is the shallow-layer low regime on Dongying Formation, and minimum speed is 1900m/s, according to formula (1), I district mesh scale is 6m; II district is small complex fault block district between Dongying Formation and substrate and hollow target area, and minimum speed is 3000m/s, and according to formula (1), II district mesh scale is 10m; III district is the mid-deep strata district under substrate and the steep major fault of height, and minimum speed is 5600m/s, and according to formula (1), III district mesh scale is 18m.

3, set man-made explosion, just drilling its propagation property in the Moving grids model of above-mentioned three kinds of different scales, determining man-made explosion wave field, obtain its border wave field, and be saved in CPU.

The man-made explosion that 3.1 settings match with three-dimensional earthquake data before superposition dominant frequency, man-made explosion wavelet is the zero phase Ricker wavelet of dominant frequency 30hz.In the Moving grids model of above-mentioned three kinds of different scales, set up the propagation property relation of man-made explosion wave field and underground medium physical parameter, the relational expression of its propagation property as shown in Equation 2.Utilize formula 2, just drilling man-made explosion from 0 moment to propagation amplitude, frequency, the phase attributes the Moving grids model being engraved in three kinds of different scales time maximum, thus determining man-made explosion wave field.In man-made explosion wavefield forward modeling process, GPU fine grained parallel calculates man-made explosion wave field propagation property, utilizes reverse-time migration big gun territory to offset incoherent feature simultaneously, and use CPU to carry out big gun territory coarse grain parallelism, common raising just drills efficiency.

$P_{i,j,k}^{t+1/2}=P_{i,j,k}^{t-1/2}+\frac{1}{\mathrm{\&rho;}}{\mathrm{\&Delta;t\&Sigma;}}_{n=1}^N\&lsqb;\frac{a_n^1}{{\mathrm{dx}}_n}{v}_{x_{i+\frac{2n-1}{2},j,k}}^t-\frac{a_n^2}{{\mathrm{dx}}_n}{v}_{x_{i-\frac{2n-1}{2},j,k}}^t\&rsqb;---\left(2\right)$

Wherein, i, j, k represent the net point coordinate in the Moving grids model of three dimensional depth territory velocity field different scale, and x, y, z representation space coordinate, t represents the time, P ^t(i, j, k) represents the normal stress wave field value at t discrete grid block (x, y, z) place, represent the speed component wave field value at t underground (x, y, z) place, represent represent Δ x, Δ y, Δ z represents the mesh spacing in (x, y, z) direction, and Δ t represents time step, and ρ represents density, and α represents difference coefficient, N representation space grid discrete number, to represent in the Difference Calculation of current the net point 2N rank forward the difference coefficient of individual grid, to represent in the Difference Calculation of current the net point 2N rank backward the difference coefficient of individual grid;

3.2 according to time order of occurrence, record the 3 D stereo wave field snapshot of the man-made explosion wave field determined by step 3.1 from 0 to each moment of maximum moment value, obtain the border wave field in the top of 3 D stereo wave field snapshot, the end, 6 faces, front, rear, left and right, and be saved in CPU.

4, hollow work area, Dong Pu oil field high-precision three-dimensional earthquake data before superposition is utilized, inverse time extrapolates its propagation property in the Moving grids model of above-mentioned three kinds of different scales, obtains the inverse time extrapolation wave field of hollow work area, Dong Pu oil field high-precision three-dimensional earthquake data before superposition.

4.1 by the wave field value assignment of single big gun maximum moment of hollow for Dong Pu oil field work area high-precision three-dimensional earthquake data before superposition in its initial extrapolation wave field;

4.2 set up three-dimensional earthquake data before superposition extrapolation wave field and the propagation property relation of underground medium physical parameter in the Moving grids model of above-mentioned three kinds of different scales, and the relational expression of its propagation property as shown in Equation 3.Utilize formula 3, initial in backward extrapolation step 4.1 successively propagation amplitude, frequency, the phase attributes of extrapolation wave field in the Moving grids model of three kinds of different scales, calculate the three-dimensional earthquake data before superposition extrapolation wave field in each moment from the maximum moment in moment to 0, and add three-dimensional earthquake data before superposition in the corresponding moment and calculate, obtain the inverse time extrapolation wave field of final three-dimensional earthquake data before superposition:

$P_{i,j,k}^{t-1/2}=P_{i,j,k}^{t+1/2}+\frac{1}{\mathrm{\&rho;}}{\mathrm{\&Delta;t\&Sigma;}}_{n=1}^N\&lsqb;\frac{a_n^1}{{\mathrm{dx}}_n}{v}_{x_{i+\frac{2n-1}{2},j,k}}^t-\frac{a_n^2}{{\mathrm{dx}}_n}{v}_{x_{i-\frac{2n-1}{2},j,k}}^t\&rsqb;---\left(3\right)$

Wherein, formula (3) is inverse time calculating on time orientation.

5, utilize step 3.2 to be kept at the border wave field of the man-made explosion wave field in CPU and the propagation property relation of described man-made explosion wave field and underground medium physical parameter, on GPU, the inverse time of the three-dimensional earthquake data before superposition that reconstruct generates with step 4 extrapolates wave field man-made explosion wave field in the same time.

6, according to zero-lag cross-correlation image-forming condition as shown in Equation 4, the inverse time extrapolation wave field of the three-dimensional earthquake data before superposition that the man-made explosion wave field reconstruct step 5 and step 4 generate carries out cross-correlation in each moment, and the cumulative cross correlation value from the maximum moment in moment to 0, determine single big gun reverse-time migration imaging results.

According to superposition sum formula as shown in Equation 5, to the cumulative summation of single big gun reverse-time migration imaging results all in hollow work area, Dong Pu oil field, and then determine the reverse-time migration imaging results of Dong Pu oil field hollow shallow-layer small complex fault block district's reverse-time migration imaging results as shown in Figure 3 and the Dong Pu oil field hollow mid-deep strata height steep fault shown in Fig. 4.

$I(x,y,z)=\frac{{\&Integral;}_0^{t_{max}}S(x,y,z,t)R(x,y,z,t)dt}{{\&Integral;}_0^{t_{max}}{S}^2(x,y,z,t)dt}---\left(4\right)$

$I_{all}(x,y,z)=\underset{i=1}{\overset{nshot}{\&Sigma;}}{I}_i(x,y,z)---\left(5\right)$

I (x, y, z) is single big gun reverse-time migration imaging results, and S is artificial source wavefield, and R is the inverse time extrapolation wave field of three-dimensional earthquake data before superposition.

I _allthe reverse-time migration imaging results in (x, y, z) work area, I _i(x, y, z) is single big gun reverse-time migration imaging results, and nshot is the total big gun number in work area.

Comparative example: utilize existing uncommon Hough pre-stack depth migration imaging technology to carry out pre-stack depth migration imaging to hollow work area, Dong Pu oil field high-precision three-dimensional earthquake data before superposition, obtains the depth migration imaging result of Dong Pu oil field hollow shallow-layer small complex fault block district depth migration imaging result as shown in Figure 5 and Dong Pu oil field hollow mid-deep strata height steep fault as shown in Figure 6.

By Fig. 3 and Fig. 5 contrast, utilize the inventive method higher to complicated minor fault imaging precision, breakpoint is more clear, features its spread in underground accurately.

By Fig. 4 and Fig. 6 contrast, utilize the inventive method better to height steep major fault cross-sectional imaging continuity, tomography both sides formation contact clearly.

Claims (6)

1. a seismic data pre-Stack Reverse formation method, is characterized in that comprising the following steps:

(1) utilize the three-dimensional earthquake data before superposition in work area to be evaluated, set up its three dimensional depth territory velocity field;

(2) utilize three dimensional depth territory, work area to be evaluated velocity field, set up the Moving grids model of three dimensional depth territory velocity field different scale;

(3) set man-made explosion, just drilling its propagation property in the Moving grids model of three dimensional depth territory velocity field different scale, determining man-made explosion wave field, obtain its border wave field, and be saved in CPU;

(4) utilize three-dimensional earthquake data before superposition, the inverse time extrapolates its propagation property in the Moving grids model of three dimensional depth territory velocity field different scale, obtains the inverse time extrapolation wave field of three-dimensional earthquake data before superposition;

(5) utilize step (3) to be kept at the border wave field of the man-made explosion wave field in CPU, on GPU, the inverse time of the three-dimensional earthquake data before superposition that reconstruct generates with step (4) extrapolates wave field man-made explosion wave field in the same time;

(6) the inverse time extrapolation wave field of the three-dimensional earthquake data before superposition that the man-made explosion wave field utilizing step (5) to reconstruct and step (4) generate, determine single big gun reverse-time migration imaging results, and then determine the reverse-time migration imaging results in work area to be evaluated.

2. a kind of seismic data pre-Stack Reverse formation method according to claim 1, is characterized in that: the Moving grids model method setting up three dimensional depth territory velocity field different scale is:

(1) physical parameter in the work area to be evaluated that the three dimensional depth territory velocity field described in utilization characterizes, carries out subregion to work area to be evaluated, determines the mesh scale of different subregion;

(2) discrete processes of the mesh scale of different subregions by determining above-mentioned steps, sets up the Moving grids model of three dimensional depth territory velocity field different scale.

3. a kind of seismic data pre-Stack Reverse formation method according to claim 1, it is characterized in that: determine man-made explosion wave field, the method obtaining its border wave field is the man-made explosion that setting matches with three-dimensional earthquake data before superposition dominant frequency, just drilling its propagation property in the Moving grids model of three dimensional depth territory velocity field different scale, determine man-made explosion wave field, obtain the border wave field of man-made explosion wave field.

4. a kind of seismic data pre-Stack Reverse formation method according to claim 2, is characterized in that: the defining method of the mesh scale of different subregion is:

$\mathrm{\&Delta;}x=\mathrm{\&Delta;}y<0.2*\frac{V_{min}}{2*feq}$

Wherein, Δ x is the discrete yardstick of grid in x direction, and Δ y is the discrete yardstick of grid in y direction, v _minfor different subregion minimum speed, feq is dominant frequency of seismic wavelet.

5. a kind of seismic data pre-Stack Reverse formation method according to claim 3, it is characterized in that: just drilling the propagation property of man-made explosion wave field in the Moving grids model of three dimensional depth territory velocity field different scale, determine man-made explosion wave field, the method obtaining the border wave field of man-made explosion wave field is:

(1) in the Moving grids model of different scale, the propagation property relation of man-made explosion wave field and underground medium physical parameter is set up, shown in the following formula of relational expression of its propagation property:

$P_{i,j,k}^{t+1/2}=P_{i,j,k}^{t-1/2}+\frac{1}{\mathrm{\&rho;}}{\mathrm{\&Delta;t\&Sigma;}}_{n=1}^N\&lsqb;\frac{a_n^1}{{\mathrm{dx}}_n}{v}_{x_{i+\frac{2n-1}{2},j,k}}^t-\frac{a_n^2}{{\mathrm{dx}}_n}{v}_{x_{i-\frac{2n-1}{2},j,k}}^t\&rsqb;$

Wherein, i, j, k represent the net point coordinate in the Moving grids model of three dimensional depth territory velocity field different scale, and x, y, z representation space coordinate, t represents the time, P ^t(i, j, k) represents the normal stress wave field value at t discrete grid block (x, y, z) place, represent the speed component wave field value at t underground (x, y, z) place, represent represent Δ x, Δ y, Δ z represents the mesh spacing in (x, y, z) direction, and Δ t represents time step, and ρ represents density, and α represents difference coefficient, N representation space grid discrete number, to represent in the Difference Calculation of current the net point 2N rank forward the difference coefficient of individual grid, to represent in the Difference Calculation of current the net point 2N rank backward the difference coefficient of individual grid;

(2) the propagation property relation of above-mentioned man-made explosion wave field and underground medium physical parameter is utilized, just drilling man-made explosion from 0 moment to propagation amplitude, frequency, the phase attributes the Moving grids model being engraved in three dimensional depth territory velocity field different scale time maximum, thus determine man-made explosion wave field, according to time order of occurrence, record the 3 D stereo wave field snapshot from the man-made explosion wave field in 0 to each moment of maximum moment value, obtain the border wave field in the top of 3 D stereo wave field snapshot, the end, 6 faces, front, rear, left and right, and be saved in CPU;

(3) in man-made explosion wavefield forward modeling process, GPU fine grained parallel calculates man-made explosion wave field propagation property, utilizes reverse-time migration big gun territory to offset incoherent feature simultaneously, uses CPU to carry out big gun territory coarse grain parallelism.

6. a kind of seismic data pre-Stack Reverse formation method according to claim 1 or 3, is characterized in that: the method setting up the inverse time extrapolation wave field of three-dimensional earthquake data before superposition is:

(1) by the wave field value assignment of single big gun maximum moment of high-precision three-dimensional earthquake data before superposition in its initial extrapolation wave field;

(2) in the Moving grids model of three dimensional depth territory velocity field different scale, set up the propagation property relation of three-dimensional earthquake data before superposition extrapolation wave field and underground medium physical parameter, the relational expression of its propagation property is as shown in following formula; Utilize following formula, propagation amplitude, frequency, the phase attributes of the backward above-mentioned extrapolation wave field initially of extrapolation in the Moving grids model of different scale successively, calculate the three-dimensional earthquake data before superposition extrapolation wave field in each moment from the maximum moment in moment to 0, and add three-dimensional earthquake data before superposition in the corresponding moment and calculate, obtain the inverse time extrapolation wave field of final three-dimensional earthquake data before superposition:

$P_{i,j,k}^{t-1/2}=P_{i,j,k}^{t+1/2}+\frac{1}{\mathrm{\&rho;}}{\mathrm{\&Delta;t\&Sigma;}}_{n=1}^N\&lsqb;\frac{a_n^1}{{\mathrm{dx}}_n}{v}_{x_{i+\frac{2n-1}{2},j,k}}^t-\frac{a_n^2}{{\mathrm{dx}}_n}{v}_{x_{i-\frac{2n-1}{2},j,k}}^t\&rsqb;$

Wherein, above-mentioned formula is inverse time calculating on time orientation.