CN101937100B - Pre-stack depth migration method - Google Patents
Pre-stack depth migration method Download PDFInfo
- Publication number
- CN101937100B CN101937100B CN201010255325A CN201010255325A CN101937100B CN 101937100 B CN101937100 B CN 101937100B CN 201010255325 A CN201010255325 A CN 201010255325A CN 201010255325 A CN201010255325 A CN 201010255325A CN 101937100 B CN101937100 B CN 101937100B
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- msup
- mover
- migration
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000013508 migration Methods 0.000 title claims abstract description 147
- 230000005012 migration Effects 0.000 title claims abstract description 147
- 238000000034 method Methods 0.000 title claims abstract description 104
- 238000003384 imaging method Methods 0.000 claims abstract description 48
- 238000004458 analytical method Methods 0.000 claims abstract description 43
- 239000013598 vector Substances 0.000 claims abstract description 22
- 230000010354 integration Effects 0.000 claims abstract description 11
- 238000006073 displacement reaction Methods 0.000 claims description 3
- 239000004576 sand Substances 0.000 claims description 2
- 238000001514 detection method Methods 0.000 claims 2
- 230000000694 effects Effects 0.000 abstract description 23
- 238000004364 calculation method Methods 0.000 abstract description 12
- 238000000926 separation method Methods 0.000 abstract description 6
- 238000010276 construction Methods 0.000 abstract 1
- 238000012545 processing Methods 0.000 description 27
- 239000010410 layer Substances 0.000 description 17
- 238000005516 engineering process Methods 0.000 description 14
- 238000004422 calculation algorithm Methods 0.000 description 13
- 238000013499 data model Methods 0.000 description 12
- 230000008569 process Effects 0.000 description 10
- 238000011160 research Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 6
- 230000006870 function Effects 0.000 description 5
- 238000004613 tight binding model Methods 0.000 description 5
- 238000004321 preservation Methods 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 230000009286 beneficial effect Effects 0.000 description 3
- 238000012937 correction Methods 0.000 description 3
- 230000006872 improvement Effects 0.000 description 3
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000009933 burial Methods 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 230000001427 coherent effect Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- 238000005755 formation reaction Methods 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 239000002245 particle Substances 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 230000002441 reversible effect Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009918 complex formation Effects 0.000 description 1
- 238000007906 compression Methods 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 230000001351 cycling effect Effects 0.000 description 1
- 238000013016 damping Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 239000011229 interlayer Substances 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003208 petroleum Substances 0.000 description 1
- 230000010363 phase shift Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 239000000523 sample Substances 0.000 description 1
- 230000035939 shock Effects 0.000 description 1
- 238000012732 spatial analysis Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000002123 temporal effect Effects 0.000 description 1
Images
Landscapes
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention relates to the field of seismic prospecting, in particular to a pre-stack depth migration method. The method comprises the following steps of: performing migration velocity analysis on seismic data; calculating travel time, arc length, emergence angle and incidence angle by ray tracing; performing migration aperture calculation; and performing pre-stack depth migration by using a Kirchhoff integration vector migration formula. The method integrates migration velocity analysis, migration aperture selection, ray tracing and Kirchhoff integration formula, does not need wave field separation, and realizes multi-component simultaneous migration and accurate homing of converted wave. Based on the practical exploration data, the pre-stack depth migration method has good surface imaging effect and high imaging resolution, the continuity of the deep reflecting in-phase axis on a converted wave imaging section is improved, and the construction is clearer.
Description
Technical Field
The invention relates to the field of seismic exploration, in particular to a prestack depth migration method.
Background
Migration is the last step in seismic data processing and is also the most important and critical process, especially for complex subsurface media.
Due to the fact that the underground medium horizon is not horizontal, the seismic stacking section cannot reflect the real horizon form, therefore, seismic data need to be subjected to migration processing on the basis of stacking, false time shift of horizon reflection generated due to the formation dip angle is corrected, reflection of the underground horizon is returned, and the real fluctuation form and the occurrence form of the formation are reflected.
The migration processing method is divided into two types, one type is carried out on the basis of superposition and is called post-stack migration; one type of migration processing is performed prior to stacking, referred to as pre-stack migration. The pre-stack migration precision is higher than the post-stack migration, and the homing capability of the complex diffracted waves is better; however, the pre-stack migration processing requires a large amount of computational effort, and has been widely developed and applied with the improvement of computer technology over the last decade. The prestack migration processing is divided into prestack time migration and prestack depth migration according to the difference of migration methods, mainly the difference of time-space domain where the methods are located. Prestack depth migration has gained greater attention in the application field because the output profile directly reflects the subsurface medium morphology of the true spatial depth domain.
There is currently very little research in the industry regarding prestack depth migration of multi-component seismic data, whether it is relatively well-established to apply post-stack time migration. The pre-stack migration processing of multi-component seismic data is an offensive and critical hotspot and difficulty in the field of seismic exploration at present, and most of current research and trial application is limited to pre-stack time migration and far from entering production practice application. Moreover, most researches are limited to the migration processing of C waves or equivalent converted waves on the basis of longitudinal wave migration along the same idea as pure longitudinal wave migration, the space vector relation among multi-component seismic signals is completely split, the migration processing of single pure waves is performed, abundant underground medium information contained in wave fields is lost, and the imaging effect is poor.
Li Ed of errors in the university of Edinburgh, UK, offered a research group (Hengchang Dai, and Xiang-Yang Li. Effect of errors in the migration coordination model of PS-converted waves on transit time access in prestack timing in week and iso regional media in geoprediction timing 2008, 73 (5): S195-S205) (Xiang-Yang Li, Jianxin Yuan. converted-wave transformed delivery-point equations in layred VTI media: Theory and seismic arrival timing in journal of Applied geoprediction timing, 2003, 34: 318) tested for successful preliminary multi-component migration timing and preliminary seismic offset component tests, which were Applied in preliminary seismic offset time experiments, but further developed for preliminary seismic component application. In addition, wearer and lie are dealing with numerical simulations and error discussions of the effect of velocity on migration results in weak anisotropic medium converted shear wave pre-stack time migration. The method respectively discusses the contributions of the converted transverse wave velocity, the vertical velocity ratio, the effective velocity ratio and the accuracy of the anisotropic parameters to the travel time calculation error, and has better reference value for improving the time domain converted wave imaging; and some ideas can be used for reference to the prestack depth migration of the converted waves.
In the last 80 th century, the concept of vector deviation was first proposed on the basis of Kuo (Kuo, J.T., and Dai, T., 1984, Kirch temporal wave deviation for the case of non-coherent source and section: Geophysics, 49(8), 1223-field 1238), and the elastic wave Kirchoff integral formula (Qinhaoh, Guo, Wang, elastic wave Kirchoff integral deviation method, geophysical report, 1988, 31 (5): 577-field 587) was developed. The method is a multi-component simultaneous migration method of elastic waves under the assumption of a uniform isotropic medium, two-component data are extended to underground imaging points in a reverse time mode through the travel time of PP waves and PS waves, and the vibration mode of three components when each imaging point is just excited is obtained. They have tested this approach to relatively good results through numerically simulated synthetic seismic records, and have initially explored it in this direction, but to date there has been no further development.
The method has the following two disadvantages: 1) the obtained vibration energy of the total two components does not realize wave field separation; 2) the offset result of this method is that the energy of a shot propagates to the shock excited for the subsurface reflection interface. The vibration vectors of different shot points when excited by the same underground reflection point are different, so that the method cannot be used for superposition of multi-shot records. Moreover, their research is an approximate prestack depth migration, i.e., the migration process is still performed in the time domain, but the time-depth conversion is performed later to convert the output migration profile into a depth profile. And for the calculation of travel time, the method only stays at the aspect of single-interface straight ray tracing, the difference of the offset aperture is not considered for different depth imaging, the offset velocity is relatively simple to obtain, the superposition velocity is simply borrowed, and then the conversion of the layer velocity is carried out, so the effect is not ideal.
The invention of the patent protection offset velocity analysis, publication number US 2009/0257308a1, is a technique that continually updates the velocity model based on the dynamic correction error so that the final velocity model approximates the true velocity of the subsurface medium. The patent is a velocity analysis technology based on a residual curvature analysis method, and the technology is based on a small offset assumption and a horizontal horizon medium assumption, namely, the residual travel can be approximated to be hyperbolic or parabolic. Under the conditions of severe speed transverse change or large offset, the algorithm based on the technology is approximately unsuitable; and the technology of the patent only processes and analyzes conventional longitudinal wave data and longitudinal wave speed.
The researches on Qi, snow, winter and the like (Qi, snow, Jiaguanhua, all megabit, Kirchhoff prestack time migration processing technology and application, geophysical prospecting calculation technology, 2009, 31 (2): 126-.
Zhang et al (Zhang, Ben, gan et al, application of integral method prestack depth migration technology in BS6 area, development of exploration geophysics, 2009, 32 (1): 48-55) introduces a mature and commercially applied prestack depth migration technology of conventional longitudinal wave seismic data and software thereof to process, test and apply seismic data of a certain area of a victory oil field, and compare the processing result with the imaging effect of the poststack migration (time) of the conventional longitudinal wave, thereby illustrating the advantages and application effect of the prestack depth migration.
Yeyueming and the like (Yeyueming, Lishuchun, all-million ambiguity, amplitude preservation prestack depth migration based on stable imaging conditions, oil geophysical exploration, 2009, 44 (1): 28-32) mainly aim at the problems of damping type deconvolution processing in the conventional longitudinal wave prestack depth migration processing, and improve the calculation stability of migration processing by Gaussian smoothing processing of the unstable calculation phenomenon caused by the denominator approaching zero, and simultaneously achieve the purpose of partial pressure noise.
The seismic data aimed at by the related technology are single-component and conventional longitudinal wave data, none of the related technologies relates to the processing of multi-component seismic data, and the related technology and theoretical method are not suitable for the processing of three-component seismic data.
The invention patent with publication number CN 101598805a only aims at the difference between the propagation time of the longitudinal wave and the converted transverse wave, and calibrates and picks up the different reflection phases of the same layer on the sections of the longitudinal wave and the converted transverse wave by the methods of velocity ratio scanning and layer tracking, so as to determine the velocity ratio of the longitudinal wave and the transverse wave of each corresponding phase, and further realize the compression of the sections of the converted transverse wave according to the propagation time of the longitudinal wave, so that the comparison between the longitudinal wave and the layer of the converted transverse wave is facilitated on the same scale of the propagation time of the longitudinal wave. The protection does not relate to imaging techniques for three-component seismic data.
The invention patent with publication number CN 1797031a aims at the problems of time consumption and low calculation efficiency of processing massive data in the prestack depth migration processing of conventional longitudinal wave data, and realizes combining multiple guns and one gun by matching and combining guns with the same phase, thereby reducing the gun set number of migration processing and improving the migration processing efficiency on the basis of not damaging the migration effect. The protection does not relate to the imaging technology of the three-component seismic data.
The invention of the patent publication No. CN 101419292A is directed to a conventional processing flow of pre-stack preprocessing, static correction, velocity analysis, dynamic correction and superposition for seismic data acquired by three components to process and image converted transverse waves. The protection content of the patent is also generally called as a conventional post-stack processing method in the field of three-component seismic data processing, and only relates to a time migration technology.
The invention patent protection of publication number US 2004/0117123a1 is directed to a technique for prestack time migration processing of conventional pure compressional data, and the technique is implemented in the frequency-wavenumber domain. The method is characterized in that seismic data are transformed into a frequency-wavenumber domain, and then linear fitting of phase shift is carried out on ray parameters of any point in a given space according to a travel time diagram, so that a prestack time migration algorithm is adapted to the complex medium imaging problem with severe transverse change of speed.
The inventive patent protection publication US 2010/0054082a1 is not related to a specific migration algorithm, more independent of multi-component prestack depth migration. In any migration process, the simulation of the wave field is involved, which is not done by ray tracing travel time calculations, but by finite difference simulations of the acoustic wave equations. The patent mainly aims at the problem that the wave field of any point of a storage space required in pre-stack reverse time migration processing causes a large amount of consumption of memory and hard disk space in migration processing, and migration efficiency is low.
Prestack depth migration is an effective tool for imaging complex formations where the dip and velocity of the earth layers vary strongly laterally. But the prestack depth migration relies on a more accurate velocity model. Conventional velocity modeling methods cannot achieve sufficient accuracy. Meanwhile, because the prestack depth migration is very sensitive to a velocity model, the prestack depth migration is used as a powerful velocity analysis tool, and two migration velocity analysis methods are developed: depth-focusing analysis (DFA) and residual-curvature analysis (RCA). Two problems must be solved to perform the offset velocity analysis: 1) establishing a standard for evaluating whether the speed is accurate or not; 2) how to perform the speed update.
The principle of depth focus analysis is to measure the velocity error from the superimposed energy. When the offset depth and the depth of focus are the same, the speed is up to the requirement. Residual curvature analysis measures the velocity error from the residual travel time difference. When the imaging depths of different shot-geophone distances are the same, the speed meets the requirement.
Because the velocity analysis is usually based on a coarse initial iteration velocity, a robust velocity analysis method can reduce the number of iterations while ensuring velocity convergence. Conventional velocity analysis methods are generally based on the following three assumptions: firstly, a layered uniform medium; secondly, small offset; and thirdly, horizontal stratum.
The DFA method of MacKay and Abma (MacKay S.and Abma R.imaging and localization with depth-focusing analysis. Geophysics.1992, 57: 1608-.
In RCA, the migration method used for velocity analysis is different, and the assumption is different: 1. shot-to-shot or geophone offset: the method proposed by Al-Yahay (Al-Yahya K.1989. spatial analysis by iterative profile analysis, Geophysics, 1989, 54: 718-; the method of Lee and Zhang (Lee W.and Zhang L.Residual shotprofile learning, Geophysics, 1992, 57: 815-. 2. Common offset: the method of Deregowski (Deregowski S.M. common-offset migration and velocity analysis. first BatTreak, 1990, 8 (6): 225-234) uses the first two assumptions.
The methods are all based on the small offset assumption, and under the condition of small offset, the rest travel time can be approximate to hyperbolic curve or parabolic curve. Using this approximation is not appropriate when the lateral variation of the speed is severe. Therefore, Liu (Liu Z. and Bleistein N. migration velocity analysis: the acoustic and iterative algorithm, Geophysics, 1995, 60: 142-.
Disclosure of Invention
The invention aims to solve the technical problem of providing a three-component prestack depth migration method based on a vector wave field aiming at multi-component seismic data, which does not need to carry out wave field separation, realizes multi-component migration and realizes accurate homing of converted waves.
The technical scheme for solving the technical problems is as follows: a pre-stack depth migration method, comprising the steps of:
step 3, calculating the offset aperture;
and 4, performing prestack depth migration by using a vector migration formula of a Kirchhoff integration method.
The invention has the beneficial effects that: the method integrates migration velocity analysis, migration aperture selection, ray tracing and Kirchhoff integral formula, does not need to carry out wave field separation, realizes multi-component migration and realizes accurate homing of converted waves.
On the basis of the technical scheme, the invention can be further improved as follows.
Further, the offset velocity analysis in step 1 adopts a co-imaging point velocity analysis method based on depth focusing.
The method has the advantages that the depth focusing-based co-imaging point speed analysis method is modified on the residual curvature analysis method, has high convergence speed and is suitable for various speed models and offset distances.
Further, the step 1 comprises the following steps:
step a1, starting a geological model, and entering step a 2;
step a2, establishing an initial speed model, and entering step a 3;
step a3, performing a pre-stack depth migration (PSDM) based on a target horizon bottom interface, and entering step a 4;
step a4, performing offset velocity analysis, and entering step a 5;
step a5, judging whether a complete speed model is established, if yes, entering step 2, and if no, entering step a 6;
step a6, carrying out PSDM based on a target horizon bottom interface, and entering step a 7;
step a7, picking up the bottom interface of the target layer, and entering step a 8;
step a8, updating the velocity model, moving the target layer to the next layer, and proceeding to step a 3.
Wherein, the complete velocity model in step a5 is: and the seismic geologic model is complete in all layers from shallow to deep and takes speed, density and layer burial depth as parameters.
Further, in step 2, a straight ray tracing method or a minimum travel time ray tracing method based on the Fermat principle is adopted.
The method has the advantages that the calculation speed is high by adopting the straight ray tracing method, and the ray blind area can be avoided by adopting the minimum travel time ray tracing method based on the Fermat principle.
Further, in the step 4, an elastic wave Kirchhoff integration method vector migration formula suitable for the non-uniform anisotropic medium is adopted to conduct prestack depth migration.
The technical scheme has the advantages that the method is suitable for an elastic wave Kirchhoff integration method vector migration formula of the non-uniform anisotropic medium, has the characteristic of wave field separation, and can separate PP waves and PS waves while migrating; the Offset result can respectively output a PP wave section and a PS wave section, which is beneficial to respectively analyzing AVO (Amplitude Versus Offset, the variation of Amplitude along with Offset) characteristics of the PP wave and the PS wave; the Green function of the elastic wave equation is used, so that the amplitude preservation of the algorithm is improved; is suitable for non-uniform anisotropic medium.
The invention provides a new offset formula based on the research of Queen Miao et al (Qin Fuhao, Guo Asia, Queen Miao, elastic wave kirchhoff integral offset method, geophysical science report 1988, 31 (5): 577-charge 587), and improves the two defects. And the PP wave and the PS wave are directly output in the migration process, the coherent enhancement effect is realized when multiple shots are superposed, and the sections of the PP wave and the PS wave are obtained.
Drawings
FIG. 1 is a first schematic diagram of a prestack depth migration method of the present invention;
FIG. 2 is a second schematic diagram of the prestack depth migration method of the present invention;
FIG. 3 is a third schematic diagram of the prestack depth migration method of the present invention;
FIG. 4 is a flow chart of a prestack depth migration method of the present invention;
FIG. 5 is a flow chart of an implementation of the offset velocity analysis of the present invention;
fig. 6A is a PP wave effect diagram after shifting a piece of 2D3C seismic data of daqing by using a Common Converted Point (Common Converted Point) gather-based post-stack time shifting method;
fig. 6B is a PS wave effect diagram after a daqing piece of 2D3C seismic data is migrated using a common CCP gather-based post-stack time migration method;
FIG. 7A is a graph of PP-wave effects after migration of a piece of 2D3C seismic data from Daqing using the method of the present invention;
FIG. 7B is a graph of the PS wave effect after migration processing of a piece of 2D3C seismic data from Daqing using the method of the present invention;
FIG. 8 is a simulated data model of the prestack depth migration method of the present invention;
FIG. 9A is a plot of a Z-component record of seismic data for the data model of FIG. 8;
FIG. 9B is a plot of a seismic data X component record for the data model of FIG. 8;
FIG. 10 is a longitudinal wave migration depth profile result obtained using the prestack depth migration method of the present invention in the data model of FIG. 8;
FIG. 11 is a converted shear wave migration depth profile result obtained using the prestack depth migration method of the present invention in the data model of FIG. 8;
FIG. 12 is a longitudinal wave migration depth profile result obtained by using acoustic wave equation migration and opening the aperture with the inclination angle in the data model of FIG. 8;
FIG. 13 is a converted shear wave migration depth profile result obtained by employing acoustic wave equation migration and opening the aperture with the inclination angle in the data model of FIG. 8;
FIG. 14 is a longitudinal wave offset depth profile result obtained by using acoustic wave equation offset and fixing the aperture in the data model of FIG. 8;
FIG. 15 is a depth profile of converted shear wave excursion obtained using the acoustic wave equation excursion and fixed aperture in the data model of FIG. 8.
Detailed Description
The principles and features of this invention are described below in conjunction with the following drawings, which are set forth by way of illustration only and are not intended to limit the scope of the invention.
Principle of method
As shown in fig. 1, the three-dimensional vector wave field continuation equation of a homogeneous isotropic medium is:
wherein
Representing the particle displacement in m/s;
represents the particle velocity in m;
denotes stress in units of N/m2(ii) a m, n-1, 2, 3 respectively represent x, y, z coordinates, in m;
is the unit radial dimension of the belt-shaped body,
the projection of a unit vector representing the ray emergence angle on the m and n coordinates respectively; t (x | x') represents the travel time of the seismic source to the underground diffraction point, and the unit is s; t is tp、tsRespectively representing the travel time of PP waves and PS waves from the diffraction point to the receiving point, wherein the unit is s; r represents the distance from the diffraction point to the receiving point and is in m; vp、VsRespectively representing the longitudinal wave speed and the transverse wave speed, and the unit is m/s; v ═ Vs/VpThe velocity ratio of the longitudinal wave and the transverse wave is shown.
Due to the free surface boundary condition, the surface stress is zero, the last two terms of equation (1) can be dropped, with
The equation is expanded, and the energy of the three components is projected to the PP component (3) and the PS component (4), respectively, as shown in FIG. 2
uPP＝u1cosθsinφ+u2sinθsinφ+u3cosφ(3)
uPS＝u1cosθcosφ+u2sinθcosφ+u3sinφ(4)
Unfolding formula (2) and then bringing the unfolded formula into formulas (3) and (4) to obtain
In two-dimensional multi-component exploration, a collapsible source is typically used, so that only P-P and P-SV waves are present on theoretical seismic records. If the vertical direction is the z-axis direction and the line measuring direction is the x-axis direction, no wave field information exists on the y-component. The equations (5) and (6) are simplified into a two-dimensional elastic wave Kirchhoff integration method vector deviation equation.
However, the actual medium is non-uniform and anisotropic, the ray paths of the PP wave and the PS wave from the underground diffraction point to the surface receiving point are no longer straight lines, and the ray paths, the travel time and the emergence angle of the PP wave and the PS wave are different. As shown in fig. 3, wherein the dotted line represents the S-wave raypath and the solid line represents the P-wave raypath. There is therefore a need for an improvement in the uniform isotropic media shift formulation. Travel time t, tp、tsRoute lp、lsAngle of emergence thetap、θsAnd equivalent longitudinal and transverse wave velocity Vp\VsRespectively obtained by ray tracing and speed analysis of the inhomogeneous anisotropic medium.
Therefore, the formulas (7) and (8) can be improved
The two equations of the above equation (9) and equation (10) are the Kirchhoff integral method vector deviation equation of the elastic wave suitable for the non-uniform anisotropic medium.
Comparing the formula of the current homogeneous isotropic Kirchhoff integral migration algorithm (11) based on the acoustic wave equation:
(9) the Kirchhoff integral migration algorithm for the vector elastic wave field in the non-uniform anisotropic medium with the formula (10) has the following advantages:
1) the method has the characteristic of wave field separation, and can separate PP and PS waves while shifting;
2) respectively outputting a PP wave section and a PS wave section by the offset result, and being beneficial to respectively analyzing AVO characteristics of the PP wave and the PS wave;
3) the Green function of the elastic wave equation is used, so that the amplitude preservation of the algorithm is improved;
4) is suitable for non-uniform anisotropic medium.
According to the principle of the method, the prestack depth migration method comprises the following implementation steps:
step 3, calculating the offset aperture;
and 4, performing prestack depth migration by using a vector migration formula of a Kirchhoff integration method.
Wherein, the offset velocity analysis in the step 1 adopts a common imaging point velocity analysis method based on depth focusing; in the step 2, a straight ray tracing method or a minimum travel time ray tracing method based on the Fermat principle is adopted; and 4, performing prestack depth migration by adopting an elastic wave Kirchhoff integration method vector migration formula suitable for the non-uniform anisotropic medium.
The flow of the implementation steps of the prestack depth migration method of the invention is shown in fig. 4.
The flow of the vector prestack depth migration algorithm is as follows:
A. the method comprises the steps of circulating all tracks of SEG-Y data (the seismic data are generally organized by taking seismic tracks as a unit and stored by adopting an SEG-Y file format, wherein the SEG-Y format is one of standard tape data formats proposed by SEG (Society of exploration geophysics Society), which is the most common format of seismic data in the petroleum exploration industry);
a) inputting a SEG-Y data, and reading the coordinates of the shot point demodulator probe through the track head (judging whether the X component or the Z component is input, wherein different travel time and Green functions are adopted for different components);
b) interpolating the rough travel time table into a grid which is the same as the offset output section;
c) cycling through all imaging points in the subsurface;
● determining whether it is within the PP wave deflection aperture;
● if the aperture is in the range, calculating the Green function of the point, reading the PP component of the imaging point superimposed by the energy corresponding to the PP wave travel time of the point;
● determining whether it is within the PS wave offset aperture;
● if the aperture is in the range, calculating the Green function of the point, reading the PS component of the energy superposed to the imaging point corresponding to the PS wave travel time of the point;
d) ending the imaging point cycle;
and B, finishing the cycle of the SEG-Y data body, writing the offset result into the file, and finishing the offset.
The general algorithm implementation flow is shown in table 1.
TABLE 1 flow chart of vector prestack depth migration algorithm
In the implementation process of the prestack depth migration method, the selection of the migration aperture and the acquisition of the migration speed are two main control factors related to the migration effect and the calculation efficiency, and the prestack depth migration method is improved as follows.
Offset aperture
The imaging quality of the Kirchhoff integration method vector migration formula for prestack depth migration can be improved by selecting a proper migration aperture. Too large an offset aperture can reduce the signal-to-noise ratio; otherwise, valid information is missed. The paths of the PP wave and the PS wave in the medium are different, so that different apertures need to be adopted.
The choice of aperture is the factor that has the greatest impact on the imaging results (preservation, accuracy). The Fresnel band offset, gaussian beam offset, finite aperture offset proposed in the prior art documents are all in fact the true position of the interfacial reflection energy is found by ray theory, opening an aperture around this point. In this process, the operations of the trigonometric function and the inverse trigonometric function need to be repeated, and a large amount of calculation is required. In addition, in the migration process, single shot data needs to be read in at one time. Therefore, in actual production, based on the assumption of horizontal interface, the aperture is opened near the CMP (Common Middle Point), and when encountering the inclination angle, the aperture is simply enlarged.
The invention adopts the same principle as the method, but improves the algorithm by using the coordinate X of the central pointmAnd the offset distance h respectively replaces the coordinates of the shot point and the receiving point, so that the calculation of a trigonometric function and an inverse trigonometric function is avoided, the operation efficiency is improved, and the method can be applied to an offset method based on an input channel.
Offset velocity analysis
The offset velocity analysis in the present invention uses a depth focusing based co-imaging point velocity analysis method based on RCA method (Liu Z. and Bleistein N. simulation analysis: Theory and an iterative algorithm, Geophysics, 1995, 60: 142-. In the common imaging point trace set, the imaging depth z is a function of the offset h. When the migration speed is equal to the underground real speed, the imaging depth is irrelevant to the offset distance, and the speed is updated on the basis of the imaging depth and the offset distance, wherein the updating formula is as follows:
wherein, δ λ is a speed updating step length, and the new iteration speed is V + δ λ; k represents a common imaging point, and M represents a shot-geophone distance; wherein,
solving by ray tracing, wherein t is interlayer travel time;
imaging depths representing different offsets; z is obtained by picking on a common imaging point gather.
Fig. 5 shows a flow of implementing offset velocity analysis, which includes the steps of:
step a1, starting a geological model, and entering step a 2;
step a2, establishing an initial speed model, and entering step a 3;
step a3, carrying out PSDM based on a target horizon bottom interface, and entering step a 4;
step a4, performing offset velocity analysis, and entering step a 5;
step a5, judging whether a complete speed model is established, if yes, entering step 2, and if no, entering step a 6;
step a6, carrying out PSDM based on a target horizon bottom interface, and entering step a 7;
step a7, picking up the bottom interface of the target layer, and entering step a 8;
step a8, updating the velocity model, moving the target layer to the next layer, and proceeding to step a 3.
Wherein, the complete velocity model in step a5 is: and the seismic geologic model is complete in all layers from shallow to deep and takes speed, density and layer burial depth as parameters.
The method has high convergence rate and is suitable for various speed models and offset distances.
FIG. 6A is a graph showing the PP wave effect after migration of a piece of Daqing 2D3C seismic data using a common CCP gather-based post-stack time migration method; FIG. 6B is a graph showing the PS wave effect after migration of a piece of Daqing 2D3C seismic data using a common CCP gather-based post-stack time migration method; FIG. 7A is a graph showing the PP-wave effect after migration of a piece of 2D3C seismic data from Daqing using the method of the present invention; fig. 7B is a graph showing the effect of PS waves after migration processing of a piece of 2D3C seismic data from daqing using the method of the present invention. As can be seen from comparison of fig. 6A, 6B, 7A, and 7B, the prestack depth migration method of the present invention has better effects than the previous method:
the imaging effect of the first and near-surface (0-500m) is good.
Second, the imaging resolution is improved. Fig. 6A and 6B show that, in the effect of the cheap PP wave and PS wave obtained by the post-stack time shift method based on the common CCP gather, only a plurality of continuous in-phase axes are visible on the time shift profile. After the prestack depth migration method is adopted, the resolution of the same phase axis is improved: the in-phase axis of the shallow layer (1000-. It can also be seen from fig. 7A and 7B that the in-phase axes are not continuous, but there are many small breaks, but this effect is not present in fig. 6A and 6B; dominant frequencies of the in-phase axis at a deep portion of 3000 m in fig. 7A and 7B are significantly increased, and a fault that is not seen on two time-shifted sections can be seen.
Thirdly, the continuity of the in-phase axis at 3000 meters on the converted wave imaging section is improved, and the structure is clearer.
FIG. 8 shows a simulated data model of the prestack depth migration method of the present invention. The model parameters are defined as: 1000m × 1000m, four layers; speed (m/s) Vp is 2500, 3000, 3500, 4000, 4500; vs 1500, 1800, 2100, 2400, 2700; density (g/cm)3) 2.250, 2.250, 2.250, 2.250, 2.250. Forward modeling multiple shot records by adopting finite difference of elastic waves; an observation system: the track spacing is 5m, 5 guns are arranged, each gun has 200 tracks, the coordinate x of an initial shot point is 100m, and the gun spacing dx is 200 m; sampling interval 0.5ms, number of sampling points: 2401. FIG. 9A is a plot of a Z-component record of seismic data for the data model of FIG. 8, and FIG. 9B is a plot of an X-component record of seismic data for the data model of FIG. 8. The migration results obtained by the prestack depth migration method of the present invention are shown in fig. 10 (longitudinal wave migration depth profile) and fig. 11 (converted transverse wave migration depth profile); FIG. 12 (compressional wave offset depth profile) and FIG. 13 (converted shear wave offset depth profile) show the results of the acoustic wave equation offset using aperture opening with tilt angle; fig. 14 (compressional wave offset depth profile) and fig. 15 (converted shear wave offset depth profile) show the results of the offset using the fixed aperture acoustic wave equation. Imaging point grid 100 x 100, dx being 10m, dz being 10 m; the offset aperture aper is 0.4, i.e. the aperture radius is the imaging point depth z times 0.4.
As can be seen from the comparison of FIGS. 10-15, the method of the present invention has the best migration effect and accurate and clear structure homing; secondly, an acoustic wave equation migration method of opening the aperture according to the inclination angle; the worst is the fixed aperture acoustic wave equation migration method.
The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents, improvements and the like that fall within the spirit and principle of the present invention are intended to be included therein.
Claims (1)
1. A pre-stack depth migration method, comprising the steps of:
step 1, performing migration velocity analysis on seismic data by adopting a depth focusing-based co-imaging point velocity analysis method, wherein the migration velocity analysis method comprises the following substeps:
step a1, starting a geological model, and entering step a 2;
step a2, establishing an initial speed model, and entering step a 3;
step a3, performing prestack depth migration based on a target layer position bottom interface, and entering step a 4;
step a4, performing offset velocity analysis, and entering step a 5;
step a5, judging whether a complete speed model is established, if yes, entering step 2, and if no, entering step a 6;
step a6, performing prestack depth migration based on the target horizon bottom interface, and entering step a 7;
step a7, picking up the bottom interface of the target layer, and entering step a 8;
step a8, when the offset speed is equal to the underground real speed, the imaging depth is independent of the offset distance, and the formula is based on the imaging depth
step 2, calculating travel time, arc length, emergence angle and incidence angle by adopting a straight ray tracing method or a minimum travel time ray tracing method based on the Fermat principle;
step 3, calculating the offset aperture;
and 4, performing prestack depth migration by using an elastic wave Kirchhoff integration method vector migration formula suitable for the non-uniform anisotropic medium, wherein the elastic wave Kirchhoff integration method vector migration formula suitable for the non-uniform anisotropic medium specifically comprises the following steps:
wherein u isPP(x, z, t) represents the Kirchhoff integral vector displacement of the longitudinal wave at the (x, z, t) point of the various heterogeneous media, uPS(x, z, t) represents the shear wave Kirchhoff integral vector displacement at the point; t = t (x | x') represents the travel time of the seismic source to the subsurface diffraction point, in units of s; t is tP，tSRespectively representing the travel time of a PP wave and a PS wave from a diffraction point to a receiving point in a unit s; vP、VSRespectively representing the longitudinal wave speed and the transverse wave speed in m/s; lP、lSRespectively showing ray paths of longitudinal waves and transverse waves; thetaP、θSRespectively showing the emergence angles of longitudinal waves and transverse waves; x' is the abscissa of the seismic source surface point;
represents the longitudinal and transverse velocities at the seismic source surface point x' at time t, V = VP/VSThe velocity ratio of the longitudinal and transverse waves is shown.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201010255325A CN101937100B (en) | 2010-08-17 | 2010-08-17 | Pre-stack depth migration method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201010255325A CN101937100B (en) | 2010-08-17 | 2010-08-17 | Pre-stack depth migration method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN101937100A CN101937100A (en) | 2011-01-05 |
| CN101937100B true CN101937100B (en) | 2012-10-03 |
Family
ID=43390528
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201010255325A Active CN101937100B (en) | 2010-08-17 | 2010-08-17 | Pre-stack depth migration method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN101937100B (en) |
Families Citing this family (35)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102156296A (en) * | 2011-04-19 | 2011-08-17 | 中国石油大学(华东) | Elastic reverse time migration imaging method by combining seismic multi-component |
| CN102914791B (en) * | 2011-08-05 | 2015-05-13 | 中国石油天然气集团公司 | Kirchhoff prestack time migration method for processing seismic data of undulating surface |
| CN102495426B (en) * | 2011-12-02 | 2014-10-22 | 北京中科联衡科技有限公司 | Kirchhoff integral seismic imaging method |
| CN102520443B (en) * | 2011-12-13 | 2014-01-29 | 北京中科联衡科技有限公司 | Method of prestack migration of diffraction waves |
| CN102590859B (en) * | 2011-12-31 | 2014-01-22 | 中国石油集团西北地质研究所 | Anisotropic reverse time migration method for quasi-P wave equation in transverse isotropy with a vertical axis of symmetry (VTI) medium |
| CN103424775B (en) * | 2012-05-18 | 2016-07-13 | 中国科学院地质与地球物理研究所 | Based near surface many focal points location positioning method that seismic wave depth migration is theoretical |
| CN103576190B (en) * | 2012-08-02 | 2016-06-08 | 中国石油天然气集团公司 | A kind of method solving reverse-time migration waveform discontinuity |
| CN102914796B (en) * | 2012-08-21 | 2013-11-13 | 北京多分量地震技术研究院 | Control method for acquiring offset speeds of longitudinal and transverse waves based on Gaussian beam |
| CN102841376A (en) * | 2012-09-06 | 2012-12-26 | 中国石油大学(华东) | Retrieval method for chromatography speed based on undulating surface |
| CN102901984B (en) * | 2012-09-29 | 2015-07-08 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | Method for constructing true earth surface dip angle trace gathers of seismic data |
| CN103777242A (en) * | 2012-10-24 | 2014-05-07 | 中国石油化工股份有限公司 | Speed discrimination method with combination of depth focusing and gather event flattening |
| US20140153365A1 (en) * | 2012-11-30 | 2014-06-05 | Chevron U.S.A. Inc. | System and method for producing local images of subsurface targets |
| CN103149592A (en) * | 2013-03-07 | 2013-06-12 | 天津城市建设学院 | Method for separating variable offset vertical seismic profile (VSP) wave fields |
| CN104375172B (en) * | 2013-08-15 | 2017-02-15 | 中国石油集团东方地球物理勘探有限责任公司 | Volcanic underlayer structural configuration correct imaging method |
| CN103698813B (en) * | 2013-12-26 | 2016-08-03 | 中国石油天然气集团公司 | A kind of ray-tracing procedure in seismic pre-stack time migration |
| WO2016008103A1 (en) * | 2014-07-15 | 2016-01-21 | 杨顺伟 | Two-dimensional turning ray integral method prestack depth migration method |
| CN104133240B (en) * | 2014-07-29 | 2017-02-01 | 中国石油天然气集团公司 | Large-scale collateral kirchhoff prestack depth migration method and device |
| CN104216016B (en) * | 2014-08-12 | 2017-01-25 | 中国石油天然气集团公司 | Converted-wave dynamic correcting method and system for longitudinal-wave constraint scanning |
| CN104216012A (en) * | 2014-08-27 | 2014-12-17 | 中国科学院地质与地球物理研究所 | Three-dimensional Born-Kirchhoff variable-step interpolation imaging method |
| CN105467439B (en) * | 2014-09-10 | 2018-08-31 | 中国石油化工股份有限公司 | Determine method of the anisotropic parameters to common imaging gather weighing factor |
| CN105891885A (en) * | 2014-10-20 | 2016-08-24 | 杨继东 | Prestack amplitude preservation focusing and imaging technology for specific geologic body underground |
| CN104459798B (en) * | 2014-11-26 | 2017-06-20 | 中国石油化工股份有限公司 | A kind of velocity modeling method based on RTM imagings |
| CN104749631B (en) * | 2015-03-11 | 2017-02-08 | 中国科学院地质与地球物理研究所 | Sparse inversion based migration velocity analysis method and device |
| CN105301646B (en) * | 2015-12-14 | 2017-08-01 | 中国科学院地质与地球物理研究所 | The method for building up and system of a kind of rate pattern integrated based on dual path |
| CN106291687A (en) * | 2016-07-21 | 2017-01-04 | 中国地质科学院地质研究所 | Anisotropy many ripples Gaussian beam pre-stack depth migration imaging method |
| CN107966734B (en) * | 2017-09-22 | 2019-04-02 | 中国地质大学（北京） | The vector denoising method of multi-component earthquake data |
| CN108710148B (en) * | 2018-05-29 | 2019-05-24 | 中国科学院地质与地球物理研究所 | The steady phase prestack depth migration method in three-dimensional dip domain and device |
| CN109085648B (en) * | 2018-07-16 | 2019-09-27 | 中国科学院地质与地球物理研究所 | Prestack depth migration method and device |
| CN109738945B (en) * | 2018-11-08 | 2021-01-19 | 成都捷科思石油天然气技术发展有限公司 | Method for directly generating construction diagram by using prestack depth migration result |
| CN109725354B (en) * | 2018-11-20 | 2020-07-10 | 中国石油天然气集团有限公司 | Anisotropic speed modeling method and system |
| CN110208853B (en) * | 2019-05-30 | 2020-07-14 | 中国地质大学（北京） | Wave equation amplitude-preserving migration method based on free interface seismic wave field derivative reconstruction |
| CN110850469A (en) * | 2019-11-20 | 2020-02-28 | 李志勇 | Imaging method for seismic channel wave depth migration based on kirchhoff product decomposition |
| CN113126152B (en) * | 2019-12-30 | 2024-06-25 | 中国石油天然气集团有限公司 | Depth domain speed model construction method and device |
| WO2024060171A1 (en) * | 2022-09-23 | 2024-03-28 | Saudi Arabian Oil Company | Method and system of imaging hydrocarbon reservoirs using adaptive aperture tapering in kirchhoff depth migration |
| CN115903043B (en) * | 2022-11-02 | 2024-03-15 | 中国矿业大学(北京) | Diffracted wave separation method and device |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101545986A (en) * | 2009-05-06 | 2009-09-30 | 匡斌 | Tridimensional integral prestack depth migration method based on maximum energy travel calculation |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20090257308A1 (en) * | 2008-04-11 | 2009-10-15 | Dimitri Bevc | Migration velocity analysis methods |
-
2010
- 2010-08-17 CN CN201010255325A patent/CN101937100B/en active Active
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101545986A (en) * | 2009-05-06 | 2009-09-30 | 匡斌 | Tridimensional integral prestack depth migration method based on maximum energy travel calculation |
Non-Patent Citations (1)
| Title |
|---|
| 宁书年等.叠前Kirchhoff积分偏移法纵横波振幅比剖面的提取与应用.《中国矿业大学学报》.2004,第33卷(第5期),495-498. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN101937100A (en) | 2011-01-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN101937100B (en) | Pre-stack depth migration method | |
| US9013956B2 (en) | Method and system for seismic imaging and earth modeling using beam tomography | |
| US6839658B2 (en) | Seismic processing with general non-hyperbolic travel-time corrections | |
| US6826484B2 (en) | 3D prestack time migration method | |
| US6002642A (en) | Seismic migration using offset checkshot data | |
| US7830747B2 (en) | Method for multi-azimuth prestack time migration for general heterogeneous, anisotropic media | |
| US8547786B2 (en) | Estimating and using slowness vector attributes in connection with a multi-component seismic gather | |
| US8174926B2 (en) | Method for wavefield separation for dual-sensor data using kirchhoff-type datuming and migration | |
| US9869783B2 (en) | Structure tensor constrained tomographic velocity analysis | |
| US11614554B2 (en) | Velocity model building for seismic data processing using PP-PS tomography with co-depthing constraint | |
| US6061301A (en) | Filtering of overburden azimuthal anisotropy effects from 3D seismic survey signals | |
| Hua et al. | Parsimonious 2D prestack Kirchhoff depth migration | |
| GB2522778A (en) | Methods and systems for optimizing generation of seismic images | |
| US10067264B2 (en) | Method of constraining seismic inversion | |
| CN107656308B (en) | A kind of common scattering point pre-stack time migration imaging method based on time depth scanning | |
| Cheng et al. | Azimuth-preserved local angle-domain prestack time migration in isotropic, vertical transversely isotropic and azimuthally anisotropic media | |
| WO2013093468A2 (en) | Full waveform inversion quality control method | |
| WO2013093467A1 (en) | Method of, and apparatus for, full waveform inversion | |
| Buland et al. | AVO inversion of Troll field data | |
| CN112462427B (en) | Multi-component seismic data amplitude-preserving angle domain common imaging point gather extraction method and system | |
| GB2503640A (en) | Quality Assurance in a Full Waveform Inversion Process | |
| GB2520124A (en) | Methods and systems for attenuating noise in seismic data | |
| CN112462428B (en) | Multi-component seismic data migration imaging method and system | |
| Luo et al. | Velocity anisotropy characteristics and pre-stack time migration imaging of the marine sedimentary strata in the South Yellow Sea Basin | |
| WANG et al. | Receiver Function Forward Modeling and Migration Based on the Wave‐Equation Finite Difference Method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20191122 Address after: 100029 Beijing city Chaoyang District Beitucheng West Road No. 19 Co-patentee after: China University of Geosciences (Beijing) Patentee after: Institute of Geology and Geophysics, Chinese Academy of Sciences Address before: 100029 Beijing city Chaoyang District Beitucheng West Road No. 19 Patentee before: Institute of Geology and Geophysics, Chinese Academy of Sciences |
|
| TR01 | Transfer of patent right |