阅读说明：本技术 一种反射结构导引的反射波层析反演方法 (Reflection wave chromatography inversion method guided by reflection structure ) 是由 王华忠 吴成梁 冯波 于 2021-10-12 设计创作，主要内容包括：本发明提供了一种基于反射结构导引的反射波层析速度反演方法,能够适用于复杂的地震数据,能够有效的筛选特征反射波数据,并基于特征反射波数据进行层析速度反演,有效地更新背景速度模型。为解决在实际数据中,由于地震数据复杂,信噪比低,导致反射波的波动理论层析速度反演不收敛问题。本发明首先基于成像剖面提取特征反射结构,然后基于特征反射结构生成观测的特征反射波和模拟的特征反射波,基于特征反射波进行走时差异测量,可以稳健地提取走时差异,更新速度模型,从而为反射波层析速度反演实用化提供一种新方法。(The invention provides a reflection wave chromatographic velocity inversion method based on reflection structure guidance, which can be applied to complex seismic data, can effectively screen characteristic reflection wave data, and can perform chromatographic velocity inversion based on the characteristic reflection wave data to effectively update a background velocity model. The method aims to solve the problem that in actual data, due to the fact that seismic data are complex and the signal-to-noise ratio is low, the inversion of the wave theoretical chromatographic velocity of reflected waves is not converged. The method comprises the steps of firstly extracting a characteristic reflection structure based on an imaging profile, then generating an observed characteristic reflection wave and a simulated characteristic reflection wave based on the characteristic reflection structure, and measuring travel time difference based on the characteristic reflection wave, so that the travel time difference can be extracted steadily, and a speed model is updated, thereby providing a new method for reflected wave chromatography speed inversion practicability.)

1. A reflection tomography inversion method guided by a reflection structure is characterized by comprising the following steps:

s1: inputting observed pre-stack seismic data;

performing offset imaging to generate an imaging gather;

leveling and correcting the imaging gather, superposing to generate a focused section, and extracting a characteristic structure;

s2: extracting a characteristic imaging gather based on the characteristic result and the imaging gather, and obtaining observed characteristic reflected wave data based on the characteristic imaging gather reverse migration;

s3: in the iteration process, shifting observed pre-stack seismic data to generate an imaging section;

extracting a characteristic structure based on the imaging section, and carrying out reverse migration to obtain simulated characteristic reflected wave data;

s4: calculating the travel time difference of the characteristic reflected waves by adopting a cross-correlation method, and constructing a corresponding adjoint source;

s5: calculating the gradient of the characteristic reflection wave chromatography inversion;

s6: carrying out post-processing on the gradient, and updating a speed model by adopting structural constraint regularization;

s7: looping the processes of S3-S6, updating the speed model, and turning to the process 8 when the termination condition is reached;

s8: the process loops S1-S7, updating the input observed characteristic reflected wave data, updating the velocity model based on the updated observed characteristic reflected wave data until a termination condition is satisfied.

2. The reflection structure guided tomographic inversion method according to claim 1, wherein: at step S1, pre-stack observed seismic data is first input, the observed seismic data being generated by the following wave equation,

d_obs(x_s,x_g,t)＝R(u(x,t))

(1)

wherein, A is a wave field equation forward operator, f (t) is a seismic source, u (x, t) is a subsurface wave field of which the wave field equation is forward, and x is { x, y, z } is a spatial point coordinate of the subsurface; v is the velocity field of the subsurface space point, R (-) is the surface sampling operator, x_s＝{s_x,s_yIs the coordinates of the shot point, x_g＝{g_x,g_yT is time, d_obsIs observed pre-stack seismic data;

and performing prestack depth migration imaging on the observed seismic data to generate an imaging gather. This process can be expressed as:

I(x,h)＝A^Td_obs (2)

where I (x, h) is the generated imaging gather, h is the extended dimension, and if an offset gather is generated, then h represents the offset of the subsurface,

if the angle domain common imaging point gather is generated, h represents the size of the underground reflection field angle, A^TThe adjoint operator represents the migration process, and taking reverse time migration as an example, the imaging process is as follows:

wherein u is_f(x, t) is a snapshot of the subsurface wavefield being performed in the initial velocity model. In the resulting imaged gather, the imaged gather is not flat due to the speed inaccuracy. We performed a flattening correction on the imaged gather, then superimposed to generate the imaged profile:

wherein, I_cor(x, h) are flattened imaging gathers; and extracting a characteristic reflection structure according to the superposition section, wherein the process is represented as:

S(x)＝f[I_cor(x)] (5)

wherein S (x) is a characteristic reflection structure, and f [. cndot. ] is a structure extraction process;

based on the extracted feature reflection structure s (x) and the imaged gather I (x, h), a feature gather is extracted, which may be represented as:

S(x,h)＝f[I(x,h)|S(x)] (6)

based on the characteristic imaging gather, carrying out reverse migration on the basis of an initial velocity model to obtain observed characteristic reflected wave data

3. The reflection structure guided tomographic inversion method according to claim 1, wherein: in step 2, the observed characteristic reflected wave calculated in step S1 is input, and in each iteration process, the simulated characteristic reflected wave is extracted, and assuming that the initial iteration number is k, the velocity model of the current iteration is v_kPerforming offset based on the current speed model, wherein the obtained offset result is I_k(x) Extracting the corresponding characteristic reflection structure as S_k(x) (ii) a Performing reverse bias based on the characteristic reflection structure, and generating simulated characteristic reflection waves by the reverse bias, wherein the simulated characteristic reflection waves are as follows:

4. the reflection structure guided tomographic inversion method according to claim 1, wherein: based on the characteristic reflected wave data extracted in the steps S1 and S2, reflected waves in the characteristic reflected waves are selected for inversion, the principle of the selected reflected waves is that the reflected waves are from primary to secondary, the reflected waves are sparse, the reflected waves corresponding to a large set of deposition layers are preferentially selected, and then the reflected waves are sequentially added into the inversion;

in reflected wave travel time tomography, it is necessary to measure the travel time difference between the observed data and the simulated data after the characteristic reflected wave pickup. The difference when taking away can be obtained through local cross-correlation, and a cross-correlation error functional is adopted, so that the method comprises the following steps:

the time lag Δ t of the reflected wave is a time shift τ by which the cross-correlation function locally reaches a maximum value, that is:

wherein, T₁,T₂Is a local time range. Having obtained the travel time residual, the satellite sources y (i, x) corresponding to the individual seismic phases_s,x_gT), expressed as:

5. the reflection structure guided tomographic inversion method according to claim 1, wherein: based on the adjoint source calculated in step S3, calculating a gradient of the characteristic reflection wave tomographic inversion, the gradient consisting of two parts, one part being a cross-correlation of the background wave field at the seismic source end and the back-transmitted adjoint primary reflection wave field, and the other part being a cross-correlation of the primary reflection wave field at the seismic source end and the back-transmitted adjoint background wave field:

wherein p is_s(x, t) and p (x, t) are respectively a background adjoint field and a primary reflection adjoint field which propagate in reverse time in a background model, u₀(x, t) and u_s(x, t) are the background wavefield and the primary reflected wavefield, respectively, at the seismic source.

6. The reflection structure guided tomographic inversion method according to claim 1, wherein:

in step S6, structure-dependent regularization is applied based on the gradient obtained in step S5, the gradient is post-processed,

the background speed model updating adopts an iterative algorithm of a gradient class to be expressed as follows:

wherein alpha is_kFor the iteration step, P is a preconditioner or regularization operator.

7. The reflection structure guided tomographic inversion method according to claim 6, wherein: and when the adopted characteristic reflection meets the conditions, repeating the process from 1 to 7, updating the input observed characteristic reflection wave data again, and continuing inversion until a high-precision speed model is obtained.

8. The reflection structure guided tomographic inversion method according to claim 1, wherein: the termination condition is that the requirement of iteration times or the requirement of error reduction is met.

9. A computer storage medium, on which a computer program is stored which, when being executed by a processor, carries out the steps of the method according to claims 1-8.

Technical Field

The invention belongs to the technical field of geophysical exploration, and particularly relates to a reflection wave chromatography inversion method guided by a reflection structure.

Background

In exploring a seismic, the full wavenumber spectrum of the subsurface model can be recovered using full waveform inversion of the seismic waves. However, the classical full waveform inversion is a strong nonlinear inversion problem, and the method is limited by various reasons and is easy to converge to a local extreme value or not converge. To reduce the non-linearity of the inversion problem, a number of variations have been proposed. The currently searched technologies mainly include the following 4 technologies:

(1) and modifying a target functional method, wherein the method is realized by modifying a conventional two-norm subtraction method and completing a more stable method. For example, a cross-correlation target functional is adopted; and carrying out inversion by methods such as normalized energy and an optimal transportation method.

(2) Differences in seismic signals are measured in the transform domain, such as in the laplace domain, or in the laplace-fourier domain, using various transform methods; inverting by adopting the travel time/phase of the seismic signal; and extracting the envelope/phase for inversion by using Hilbert transform.

(3) By introducing different regularization or prior information, the stability of convergence is improved, and a result more in line with geological logic is obtained, for example, Gihonov regularization is adopted; l1 regularization; TV regularization; and sparse constraints for various models, etc.

(4) Performing inversion on the data by dividing the data and frequency bands, for example, performing inversion by adopting a strategy from low frequency to high frequency; dividing the transmitted wave and reflected wave data for inversion, and the like; performing inversion by adopting a skeleton function; and performing inversion by using the seismic data subset.

At present, the transmitted wave component of shot concentration is mainly used for updating the background speed. However, transmitted waves require large offset seismic data and the depth range of velocity updates is shallow. One more efficient method is to update the mid-depth background velocity information with the reflected wave data. Considering the tomographic velocity inversion using reflected waves, it is necessary to remove high wave number information in the gradient and generate simulated reflected wave data using an anti-migration operator. Then, in the actual seismic data, the signal-to-noise ratio is low, and the coherent and incoherent noises are more contained, so that the reflection wave chromatography inversion is not converged and cannot be applied in the actual application process. Aiming at the practicability of the reflection wave chromatography inversion, the method adopts a reflection wave chromatography method based on reflection structure guidance, avoids the problem of cycle jump by reasonably selecting characteristic data, and improves the quality of the gradient.

Disclosure of Invention

The invention aims to provide a reflection wave chromatographic velocity inversion method based on reflection structure guidance, which can be applied to complex seismic data, can effectively screen characteristic reflection wave data, and can perform chromatographic velocity inversion based on the characteristic reflection wave data to effectively update a background velocity model. The method aims to solve the problem that in actual data, due to the fact that seismic data are complex and the signal-to-noise ratio is low, the inversion of the wave theoretical chromatographic velocity of reflected waves is not converged. Firstly, a characteristic reflection structure is extracted based on an imaging profile, then an observed characteristic reflection wave and a simulated characteristic reflection wave are generated based on the characteristic reflection structure, travel time difference measurement is carried out based on the characteristic reflection wave, travel time difference can be extracted steadily, and a speed model is updated, so that a new method is provided for reflected wave chromatography speed inversion practicality.

In a first aspect,

the invention provides a reflection wave chromatography inversion method guided by a reflection structure, which comprises the following steps:

s1 inputting observed pre-stack seismic data;

performing offset imaging to generate an imaging gather;

leveling and correcting the imaging gather, superposing to generate a focused section, and extracting a characteristic structure;

s2, extracting a characteristic imaging gather based on the characteristic result and the imaging gather, and reversely biasing based on the characteristic imaging gather

Moving to obtain observed characteristic reflected wave data;

s3, shifting the observed pre-stack seismic data in an iteration process to generate an imaging section;

extracting a characteristic structure based on the imaging section, and carrying out reverse migration to obtain simulated characteristic reflected wave data;

s4: calculating the travel time difference of the characteristic reflected waves by adopting a cross-correlation method, and constructing a corresponding adjoint source;

s5: calculating the gradient of the characteristic reflection wave chromatography inversion;

s6: carrying out post-processing on the gradient, and updating a speed model by adopting structural constraint regularization;

s7: looping the processes of S3-S6, updating the speed model, and turning to the process 8 when the termination condition is reached;

s8: the process loops S1-S7, updating the input observed characteristic reflected wave data, updating the velocity model based on the updated observed characteristic reflected wave data until a termination condition is satisfied.

In the above method, the present disclosure is further limited:

the pre-stack observed seismic data is first input, which may be considered to be generated by the following wave equation,

d_obs(x_s,x_g,t)＝R(u(x,t))

(18)

wherein, A is a wave field equation forward operator, f (t) is a seismic source, u (x, t) is a subsurface wave field of which the wave field equation is forward, and x is { x, y, z } is a spatial point coordinate of the subsurface; v is the velocity field of the subsurface space point, R (-) is the surface sampling operator, x_s＝{s_x,s_yIs the coordinates of the shot point, x_g＝{g_x,g_yT is time, d_obsIs observed pre-stack seismic data;

and performing prestack depth migration imaging on the observed seismic data to generate an imaging gather. This process can be expressed as:

I(x,h)＝A^Td_obs (19)

where I (x, h) is the generated imaging gather, h is the extended dimension, and if an offset gather is generated, then h is

Representing the underground offset, if the angle domain common imaging point gather is generated, h represents the size of the underground reflection field angle, A^TIs an adjoint operator, which represents the migration process, and taking reverse time migration as an example, the imaging process can be represented as:

wherein u is_f(x, t) is a snapshot of the subsurface wavefield being performed in the initial velocity model. In the generated image trace set, due toThe speed is inaccurate resulting in an uneven gather of imaged traces. We performed a flattening correction on the imaged gather, then superimposed to generate the imaged profile:

wherein, I_cor(x, h) are flattened imaging gathers; and extracting a characteristic reflection structure according to the superposition section, wherein the process is represented as:

S(x)＝f[I_cor(x)] (22)

wherein S (x) is a characteristic reflection structure, and f [. cndot. ] is a structure extraction process;

based on the extracted feature reflection structure s (x) and the imaged gather I (x, h), a feature gather is extracted, which may be represented as:

S(x,h)＝f[I(x,h)|S(x)] (23)

based on the characteristic imaging gather, carrying out reverse migration on the basis of an initial speed model to obtain observed characteristics

Reflected wave data

In step 2, the observed characteristic reflected wave calculated in step S1 is input, and in each iteration process, the simulated characteristic reflected wave is extracted, and assuming that the initial iteration number is k, the velocity model of the current iteration is v_kPerforming offset based on the current speed model, wherein the obtained offset result is I_k(x) Extracting the corresponding characteristic reflection structure as S_k(x) (ii) a Performing reverse bias based on the characteristic reflection structure, and generating simulated characteristic reflection waves by the reverse bias, wherein the simulated characteristic reflection waves are as follows:

based on the characteristic reflected wave data extracted in the steps S1 and S2, reflected waves in the characteristic reflected waves are selected for inversion, the principle of the selected reflected waves is that the reflected waves are from primary to secondary, the reflected waves are sparse, the reflected waves corresponding to a large set of deposition layers are preferentially selected, and then the reflected waves are sequentially added into the inversion;

in reflected wave travel time tomography, it is necessary to measure the travel time difference between the observed data and the simulated data after the characteristic reflected wave pickup. The difference when taking away can be obtained through local cross-correlation, and a cross-correlation error functional is adopted, so that the method comprises the following steps:

the time lag Δ t of the reflected wave is a time shift τ by which the cross-correlation function locally reaches a maximum value, that is:

wherein, T₁,T₂Is a local time range. Having obtained the travel time residual, the satellite sources y (i, x) corresponding to the individual seismic phases_s,x_gT), expressed as:

based on the adjoint source calculated in step S3, calculating a gradient of the characteristic reflection wave tomographic inversion, the gradient consisting of two parts, one part being a cross-correlation of the background wave field at the seismic source end and the back-transmitted adjoint primary reflection wave field, and the other part being a cross-correlation of the primary reflection wave field at the seismic source end and the back-transmitted adjoint background wave field:

wherein p is_s(x, t) and p (x, t) are respectively a background adjoint field and a primary reflection adjoint field which propagate in reverse time in a background model, u₀(x, t) and u_s(x, t) are the background wavefield and the primary reflected wavefield, respectively, at the seismic source.

In step S6, structure-dependent regularization is applied based on the gradient obtained in step S5, the gradient is post-processed,

the background speed model updating adopts an iterative algorithm of a gradient class to be expressed as follows:

wherein alpha is_kFor the iteration step, P is a preconditioner or regularization operator.

And when the adopted characteristic reflection meets the conditions, repeating the process from 1 to 7, updating the input observed characteristic reflection wave data again, and continuing inversion until a high-precision speed model is obtained.

The termination condition is that the requirement of iteration times or the requirement of error reduction is met.

In a second aspect, the present invention also provides a computer storage medium having a computer program stored thereon, which when executed by a processor, performs the steps of the method of the first aspect.

Compared with the prior art, the method has the beneficial effects that

Firstly, feature reflected wave data of observation is screened through migration-reverse migration of an imaging trace gather domain, a feature imaging trace gather is extracted, then reverse migration is carried out on the basis of the feature imaging trace gather, expected feature reflected wave data can be obtained, travel time information which is the same as that of original observed seismic records can be obtained, and therefore the method can be used for subsequent velocity updating.

In addition, through the migration and the reverse migration, characteristic observation data are constructed, and the influences of noise, wavelet inconsistency and observation system irregularity can be eliminated.

Secondly, extracting characteristic structure information, screening reflection axes, and selecting flexibility of accompanying sources; and in the inversion process, effective information is added in sequence to realize stable speed updating.

Thirdly, the method has no cycle jump problem and has low dependence on the initial model. The gradient quality is high, and the method can be suitable for complex conditions.

And fourthly, structural constraint is added in the speed updating, so that the inverted speed model has geological significance and a more reliable result is obtained.

Fifth, the method has a very wide application value. The method has the advantages that the problems of low signal-to-noise ratio, strong noise, difficult identification of effective reflection axes and the like exist in marine actual data and onshore actual data, and a better inversion result can be obtained by the method aiming at the data, so that the method has wide application prospect and can provide a new method for reflected wave chromatography speed inversion practicability.

Drawings

FIG. 1 is a flow chart of a reflection wave tomographic velocity inversion based on reflection structure guidance.

Fig. 2 is a true velocity model.

FIG. 3a is observed characteristic reflected wave data extracted in sequence;

FIG. 3b is simulated characteristic reflected wave data extracted in sequence;

figure 3c is the corresponding satellite source.

FIG. 4a is a graph of velocity update amount;

fig. 4b is the velocity result of the inversion.

Fig. 5 is an error fit curve.

Figure 6 is a comparison of a longitudinally drawn line,

the dotted line is the initial velocity, the dotted line is the real velocity, and the black solid line is the inversion result.

FIG. 7 is shot gather comparison results;

wherein, the left shot gather is the observed seismic record, and the middle shot gather is the seismic record adopting initial velocity simulation; the seismic records using inverted velocity simulations are on the right.

FIG. 8 is a shot gather pull comparison result for different iterations; the first line is the correct reflected wave position, and the results of different iteration times are displayed in sequence.

Fig. 9a is an accurate velocity model.

FIG. 9b is the velocity results of the inversion.

FIG. 10a is a diagram of an extracted feature imaging gather.

FIG. 10b is an observed seismic record;

FIG. 10c shows characteristic reflections of the selected observations.

FIG. 11a is an offset imaging under a background velocity model.

Figure 11b is the corrected offset imaging result,

figure 11c shows the corresponding feature of the initial velocity model,

fig. 11d is a simulated characteristic reflected wave based on imaging profile reverse offset extraction.

FIG. 12a is a comparison of vertical line-drawing results; the dotted line is the initial velocity, the dotted line is the true velocity, and the black solid line is the inversion result.

FIG. 12b shows the results of the cross-hatch comparison. The dotted line is the initial velocity, the dotted line is the true velocity, and the black solid line is the inversion result.

Figure 13a is a graph of imaging results using an initial velocity model offset,

FIG. 13 b: shifting an imaging result by adopting an accurate speed model;

FIG. 13 c: and (5) adopting an inversion speed migration imaging result.

FIG. 14a is a plot of offset gather results from using initial velocity offset;

FIG. 14b is the offset gather results using inverted velocity offsets.

Detailed Description

The following description of specific embodiments of the present invention is provided in connection with examples to facilitate a better understanding of the present invention. Better understanding of

Example one

The reflection wave time-lapse tomography inversion effect of the invention is illustrated by a four-layer layered model. The true velocity model is shown in fig. 2. A constant speed of 1800m/s was used as the initial speed. Four axes in the imaging profile are selected as the characteristic reflection structures, so the characteristic reflection waves that generate the corresponding observations are also 4 reflection axes. In the iterative process, 4 reflection wave in-phase axes are selected to participate in the inversion at the same time. The travel time difference is measured by a cross-correlation method, wherein fig. 3a is a sequence of extracted observed characteristic reflected waves, fig. 3b is a sequence of extracted simulated characteristic reflected waves, and the accompanying source obtained based on the calculated travel time difference is shown in fig. 3 c. Since the characteristic reflected waves are sequentially extracted in the experiment, the method does not have the condition of cycle skip. The velocity update amount and inversion results for 49 iterations are shown in fig. 4a and 4 b. The corresponding convergence curve is shown in fig. 5. The results of the velocity inversion are shown in figure 6 for the vertical line-drawing. It can be seen that the velocity model is updated effectively, and other reflective layers are updated better, except that there is no effective update (no reflection axis) below the reflective interface of the deepest layer. FIG. 7 shows the result of a shot-gather comparison. The left shot gather is an observed seismic record, and the middle shot gather is a seismic record adopting initial velocity simulation; the seismic records using inverted velocity simulations are on the right. Comparing the shot gather results, it can be seen that the reflected wave is basically consistent with the real reflected wave position when the iteration is carried out for the last time. The result of the line drawing of the shot gather is shown in FIG. 8. It can also be seen that the results of the inversion produce a position of the axis of reflection that is very close to the true axis of reflection.

Example two

And further illustrating the effect of the characteristic reflected wave time-lapse tomography inversion by adopting a more complex whale type model. The true velocity model is shown in fig. 9 a. The initial velocity model was taken as a constant velocity of 2000 m/s. Selecting 3 reflection wave event axes to participate in inversion (not including the middle salt dome part), wherein fig. 10a is an extracted characteristic imaging gather, and after the reverse migration of an imaging gather domain, the result is as shown in fig. 10c, fig. 10b is an observed seismic record, comparing fig. 10b with fig. 10c, after screening, only main characteristic reflection wave data are reserved, and the travel time of the characteristic reflection wave data is consistent with that of the original observed seismic record, which shows that the data can be better screened by adopting the reverse migration of the imaging gather domain guided by the structure. In the iterative process, the simulated characteristic reflected wave data is shown in fig. 11. Wherein fig. 11a is the offset imaging result at the background speed, the offset result is not focused due to the inaccurate background speed, the focused imaging result is shown in fig. 11b after the correction, the corresponding extracted feature structure is shown in fig. 11c, and the anti-offset result of the imaging domain based on the feature structure is shown in fig. 11 d. After the inversion iteration, the velocity results of the inversion are shown in fig. 9 b. The results of the comparison of the longitudinal and transverse direction sketches are shown in figure 12. It can be seen that after the characteristic reflection tomography inversion based on the characteristic structure guidance, the velocity model is substantially close to the real velocity model. The initial velocity model migration imaging results are shown in fig. 13a, the accurate velocity model migration imaging results are shown in fig. 13b, and the inversion velocity migration imaging results are shown in fig. 13 c. The comparison shows that in the migration result obtained by the initial velocity, the position of the reflection axis is incorrect, the imaging result is unfocused, and in the result obtained by the inversion velocity, the position of the reflection axis is corrected and basically consistent with the migration result obtained by the accurate velocity model. The corresponding imaging gather results are shown in figure 14. And (3) adopting an initial velocity migration result, causing the imaging gather to bend upwards due to too low migration velocity, and flattening the imaging gather after tomographic inversion. The accuracy of the offset velocity obtained by inversion is illustrated, and the effectiveness of the chromatographic inversion is proved.

While the foregoing is directed to the preferred embodiment of the present invention, it will be understood by those skilled in the art that various changes and modifications may be made without departing from the spirit and scope of the invention.