Three-dimensional kirchhoff integral method prestack time migration fast imaging method

Abstract

The application provides a three-dimensional Kirchhoff integral method prestack time migration fast imaging method, comprising the following steps: step 1, collecting a velocity model and seismic data for migration imaging; step 2, calculating travel time of different layers of scattering points on different threads; step 3, based on an instruction set, calculating travel of four longitudinal sampling points of each trace in parallel and storing on an SSE register; step 4, for each scattering point, performing geometric diffusion compensation of amplitude energy, completing scanning and superposition of diffraction energy, and realizing migration imaging processing based on the scattering point; and step 5, performing migration imaging in a migration aperture range and completing fast imaging processing of multiple threads in a single machine. The three-dimensional Kirchhoff integral method prestack time migration fast imaging method provides a more rapid and attractive prestack time migration technology, and provides a basis and guarantee for subsequent field data acquisition site quality monitoring.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of exploration geophysics, in particular to a three-dimensional Kirchhoff integral pre-stack time migration fast imaging method. BACKGROUND

[0002] With the continuous deepening of oil and gas exploration and development, its difficulty and complexity are increasing. Seismic exploration technology needs to be developed for geological targets, and the targets are becoming more and more complex, so higher requirements are put forward for the signal-to-noise ratio, resolution, fidelity and imaging accuracy of seismic exploration data. The current conventional seismic technology has been difficult to meet these geological and engineering needs, and it is urgent to improve the accuracy of seismic exploration and the ability to solve complex geological problems, and the primary task is to obtain high-precision and high-fidelity seismic exploration original data source.

[0003] In order to improve the exploration accuracy of lithologic oil and gas reservoirs and complex targets, high-precision and high-density seismic exploration methods have been developed in recent years at home and abroad. Through the form of single-point non-combination, field data acquisition is realized. The characteristics of this exploration method are small bin, super multi-channel and large dynamic range. Compared with conventional seismic exploration, it has the advantages of high spatial sampling rate, improved frequency band, imaging accuracy and resolution. Another problem brought by this acquisition method is that the amount of data collected is huge, and a piece of three-dimensional data is measured in TB as the basic unit, and the data of field construction per day also has hundreds of GB size. How to monitor the quality of massive seismic data on site to ensure the quality of collected data has become a research hotspot and difficulty.

[0004] In the past, due to the immaturity of pre-stack migration theory and the limitation of computer capacity, post-stack migration has always dominated in seismic data processing. However, when there are steeply dipping reflectors underground, the horizontal stacking theory becomes unreliable, and CMP stacking is not common reflection point stacking. Only the structure of CIP (common imaging point) stack stacking is the best imaging to determine the imaging velocity and reflection interface position simultaneously. Therefore, pre-stack migration method is the first choice for complex structure imaging. The pre-stack migration theory cancels the assumption of input data bit zero shot, avoids the distortion caused by NMO correction after stacking, and preserves more pre-stack seismic information than post-stack migration. Compared with post-stack migration, pre-stack time migration can transfer the reflection wave energy existing in each recording channel to its real underground position under the condition that the lateral velocity variation is not very severe. Therefore, under the condition that the strata dip angle is large and the lateral velocity variation is not very severe, pre-stack time migration can achieve good imaging effect. In areas where the lateral velocity variation is not very significant and the velocity does not have sudden change, pre-stack time migration can be used instead of pre-stack depth migration to realize the migration imaging processing of data.

[0005] Kirchhoff integral migration method and micro-classification method based on wave field depth recursive extrapolation, namely the so-called wave equation method. Relative to the wave imaging method, Kirchhoff integral migration method is based on ray theory, and the scattered wave field and its normal derivative on the observation surface are integrated to realize the wave field back propagation. The geometric meaning of Kirchhoff migration is more intuitive, and the algorithm implementation adopts the diffraction superposition method, which is simple and efficient, so it is still widely used in the enterprise. Considering the timeliness of mass seismic data processing, a new three-dimensional Kirchhoff prestack time migration fast imaging technology is developed under the premise of ensuring imaging accuracy. The imaging process is optimized by travel time fast calculation of coarse and fine grid combination, multi-thread parallel and instruction set parallel, etc. On the basis of ensuring imaging effect, only relying on a single workstation, the three-dimensional mass seismic data is realized on-site fast imaging, and the purpose of quickly monitoring the quality of field acquisition data by processing means is achieved.

[0006] In the Chinese patent application with application number: CN201110036627.9, a high-precision prestack domain least square migration seismic imaging technology is involved. The steps are as follows: an initial velocity model is given, the travel time table is calculated according to the ray tracing method and stored in the memory for subsequent migration process; at the same time, the forward wave field is calculated by using high-frequency approximation theory; then the residual error between the forward wave field and the recorded wave field is calculated, and when the residual error is less than the given error standard, the iteration is stopped and the imaging result is output; when the deviation is greater than the given error standard, the difference data is migrated and the correction coefficient is calculated, and finally the imaging result is updated, and then the above operation is repeated until the calculation error is less than the given error standard, and the imaging result is output.

[0007] In the Chinese patent application with application number: CN201510424457.X, a multi-component joint Gaussian beam prestack reverse time migration imaging method is involved, which includes: inputting an initial P-S wave velocity field; reading in P-wave and converted wave seismic shot record, and determining parameters such as frequency band width, beam width, beam center interval, etc.; constructing the underground elastic vector wave field; imaging by using cross-correlation imaging condition; migrating and stacking all shot records, and finally obtaining the elastic wave imaging result.

[0008] In the Chinese patent application with the application number CN201010255325.6, a pre-stack depth migration method is disclosed. The method comprises: performing migration velocity analysis on seismic data; calculating travel time, arc length, exit angle and incident angle by ray tracing; performing migration aperture calculation; and performing pre-stack depth migration by using Kirchhoff integral method vector migration formula. The method combines migration velocity analysis, migration aperture selection, ray tracing and Kirchhoff integral formula, and does not need to perform wave field separation, and realizes multi-component simultaneous migration and accurate positioning of converted waves.

[0009] The prior art above has great difference from the present application, and cannot solve the technical problems we want to solve. Therefore, we invent a new three-dimensional Kirchhoff integral method pre-stack time migration fast imaging method. SUMMARY

[0010] The present application aims to provide a three-dimensional Kirchhoff integral method pre-stack time migration fast imaging method which solves the problem that only a single processing workstation in the field cannot perform three-dimensional fast pre-stack migration.

[0011] The object of the present application can be achieved by the following technical measures: a three-dimensional Kirchhoff integral method pre-stack time migration fast imaging method, which comprises:

[0012] Step 1: collecting a velocity model and seismic data for migration imaging;

[0013] Step 2: completing the calculation of travel time of different layers of scattering points on different threads;

[0014] Step 3: based on the instruction set, the calculation of travel of 4 longitudinal sampling points of each trace is realized in parallel, and is stored on the SSE register;

[0015] Step 4: for each scattering point, the geometric diffusion compensation of amplitude energy is performed, the scanning and superposition of diffraction energy are completed, and the migration imaging processing based on the scattering point is realized;

[0016] Step 5: performing migration imaging in the migration aperture range, and completing the fast imaging processing of multi-thread in a single machine.

[0017] The object of the present application can also be achieved by the following technical measures:

[0018] In step 1, the seismic data is processed, including noise suppression, multiple wave attenuation and amplitude energy compensation, so as to prepare a good input trace set for the next fast imaging, and the trace set is a common shot trace set or a common center point trace set.

[0019] In step 1, a velocity model required for the imaging process is constructed. This model uses the stacking velocity obtained after velocity analysis as the initial migration velocity, and after migration velocity analysis, a migration velocity model for time migration is obtained. Seismic gathers are then evenly distributed to different CPU threads to complete the data preparation work for fast time migration.

[0020] In step 2, during the Kirchhoff integral method pre-stack rapid imaging process, the subsurface strata are considered to consist of many scattering points. For each scattering point, the travel time is calculated using a double square root equation. During migration, it is assumed that the wave propagates in a straight line from the source point to the scattering point. The total travel time t is the travel time t from the source point to the diffraction point. s and the travel time t from the diffraction point to the receiving point r The sum, the formula is as follows:

[0021] t = t s +t r

[0022] First, assuming the velocity V is constant, expanding the above equation yields the double square root equation:

[0023]

[0024] Where: z0 is the depth of the scattering point, x is the position of the midpoint of the shot-receiver distance MP relative to the scattering point SP located at x=0, h is half of the offset distance from the source point to the receiver point, and v is the velocity;

[0025] To approximate the case where the velocity varies strongly with depth but weakly with lateral movement, the double square root equation can be modified as follows:

[0026]

[0027] Where the offset velocity v mig It is the approximate root-mean-square velocity estimated by Taner and Koehler at t0, where time t0 = t (x = 0, h = 0) is expressed using the average velocity v. ave Calculated two-way zero-shot-receiver distance time:

[0028]

[0029] In step 2, when calculating the time offset, according to a certain sampling interval Δt, the coarse imaging grid points are first divided, and the seismic travel time is calculated on the coarse grid using the double square root formula. Then, the travel time of the intermediate fine grid points is obtained by using linear interpolation or cubic convolution interpolation.

[0030] In step 2, the selection of △t is related to the main frequency f0 of the seismic wave, and under the premise that the longitudinal velocity does not change dramatically, the following formula is used to select the △t interval, which can basically ensure a good balance between the calculation efficiency and the imaging accuracy:

[0031] 2 / f0≤△t≤4 / f0

[0032] In the interval of two △t, the travel time of the sample points is calculated by using the double square root equation, and in the interval, the travel time is calculated by the idea of sample point interpolation. Through the combination of the two, the calculation efficiency of the travel time is improved.

[0033] In step 3, under the condition of a single workstation on site, the calculation and storage are completed by relying on the acceleration ability of the CPU, the Single Instruction Multiple Data Stream (SSE) technology is introduced, and the processor with Intel SSE instruction set support has 8 128-bit registers. Each register can store 4 32-bit single-precision floating-point numbers, and these numbers can be subjected to arithmetic and logical operations in these registers.

[0034] In step 3, after using the SSE technology, the algorithm for calculating the travel time based on the double square root is as follows:

[0035] For every 4 elements in the array, the 4 elements refer to the travel time of the scattering points at different longitudinal sampling depths. The 4 longitudinal sampling points in the array are loaded into an 128-bit SSE register; the double square root travel time of the 4 sampling numbers is calculated in one CPU instruction execution period; and the obtained 4 results are taken out and written into the memory.

[0036] Without considering the overhead of loading floating-point data in the SSE register, the algorithm calculates the travel time of 4 longitudinal sampling points at a time in one CPU period, so the calculation efficiency can be improved by about 4 times. Through the combination of instruction set parallelism, the on-site rapid imaging capability of a single workstation is improved.

[0037] In step 4, for each scattering point, the geometric diffusion compensation of the amplitude energy is performed, the scanning and stacking of the diffraction energy are completed, and the scattering point-based migration imaging processing is realized. Based on the weighted function formula of three-dimensional Kirchhoff pre-stack time migration, the compensation of the geometric diffusion loss energy of the scattering point is realized, and the amplitude preservation of the rapid imaging algorithm is improved.

[0038] In step 4, for two-dimensional,

[0039]

[0040] For three-dimensional,

[0041]

[0042] where Q m represents a weight function, z represents the depth at the imaging point, v mig represents the migration velocity, t shot represents the travel time of the source to the scattering point, t rec represents the travel time of the receiver to the scattering point.

[0043] In step 5, in the three-dimensional case, based on the conical migration aperture, the diffraction hyperboloid energy stack of each scattering point ψ xy is completed, and the migration homing processing of each scattering point is realized; the migration aperture is a circular cone with a vertex above the ground surface, and the definition is as follows: ψ xy represents the central position of the scattering point, and θ represents the vertex angle of the circular cone; the aperture is calculated according to:

[0044]

[0045] is realized, wherein v mig represents the migration velocity, and t0 represents the two-way travel time;

[0046] The imaging results of n threads are combined to complete the fast imaging processing of the entire data body.

[0047] The three-dimensional Kirchhoff integral pre-stack time migration fast imaging method in the application relates to how to realize fast imaging processing of massive three-dimensional data collected on site based on the Kirchhoff integral pre-stack time migration fast imaging technology under the limitation of only a single on-site processing workstation in the field construction site, and realize a fast judgment of the quality of the on-site collected data through the imaging results, so as to achieve the purpose of on-site quality monitoring of the on-site collected data.

[0048] Under the condition that there is no parallel machine and large processing equipment in the current field construction site, the application develops an efficient and fast imaging algorithm, only uses a single on-site processing workstation, completes the fast imaging processing of three-dimensional massive data, realizes the evaluation of the quality of the on-site collected data, and improves the quality monitoring precision of the on-site field.

[0049] The conventional Kirchhoff integral time migration algorithm is improved, under the condition limitation of a single workstation, the imaging algorithm is improved, based on the idea of migration of each input trace, the limitation of the input data trace set is broken through, the calculation of coarse and fine grid fast travel time, multi-thread parallel and amplitude seeking of the amplitude-preserving weight function and other strategies are used, the relative amplitude fidelity processing is realized, the calculation efficiency of the imaging algorithm is improved, and the pre-stack fast time migration imaging of massive seismic data is realized.

[0050] 1. The input data of the invented fast imaging algorithm is single-channel offset calculation, which has no constraint and limitation for the gather, and can be single channel of shot gather or CMP gather, etc. Based on the single-channel offset imaging idea, the imaging algorithm is free from the gather limitation, and the imaging input is more flexible.

[0051] 2. The travel time is calculated based on the double square root equation on the coarse grid, and the travel time is calculated by interpolation on the fine grid, so as to realize the calculation of fast imaging travel time and improve the calculation efficiency of the whole offset process.

[0052] 3. Based on the multi-thread parallel idea, the calculation efficiency of the algorithm is further improved. Since Kirchhoff migration is a one-to-many data mapping process, offset is performed on each imaging line of the three-dimensional imaging space, and only one channel is calculated each time, which well solves the problem that the insufficient memory overhead of multi-thread parallelism has adverse effect on the calculation efficiency. The input seismic data is shared between different threads, and the calculation performance of multi-core CPU in pre-stack time migration process is effectively improved.

[0053] 4. Through the acceleration strategy of SSE instruction set parallel, the travel time calculation and storage of four sampling points are completed at one CPU instruction execution cycle, so as to further improve the calculation efficiency of the algorithm, and finally realize the overall fast migration imaging of data.

[0054] The method is complete and scientific, and is an excellent pre-stack time migration fast imaging method, which solves the problem that three-dimensional fast pre-stack migration cannot be performed on a single processing workstation in the field, provides a more rapid and attractive pre-stack time migration technology, and provides a basis and guarantee for subsequent field quality monitoring of collected data. BRIEF DESCRIPTION OF DRAWINGS

[0055] Figure 1 The flowchart of a specific embodiment of the invented three-dimensional Kirchhoff integral method pre-stack time migration fast imaging method;

[0056] Figure 2 The path schematic diagram of Kirchhoff pre-stack time migration in a specific embodiment of the invention;

[0057] Figure 3 The travel time calculation schematic diagram of coarse and fine grid combination in a specific embodiment of the present invention;

[0058] Figure 4 The schematic diagram of three-dimensional Kirchhoff integral method migration aperture in a specific embodiment of the invention;

[0059] Figure 5 The schematic diagram of inline 305 line velocity and CMP gather in a specific embodiment of the invention;

[0060] Figure 6 This is a schematic diagram comparing the effects of different offset apertures in a specific embodiment of the present invention;

[0061] Figure 7 This is a schematic diagram of the inline 305 line pre-stack fast imaging processing profile in a specific embodiment of the present invention;

[0062] Figure 8 This is a schematic diagram of the inline305 line commercial software processing profile in a specific embodiment of the present invention;

[0063] Figure 9 This is a schematic diagram comparing the imaging effects of different volume offset data along the inline305 line in a specific embodiment of the present invention.

[0064] Figure 10 This is a schematic diagram of the CMP gather and offset velocity field of an inline 350 line in a specific embodiment of the present invention;

[0065] Figure 11 This is a schematic diagram of an inline 350-line pre-stack rapid imaging profile and a commercial software processing profile in a specific embodiment of the present invention.

[0066] Figure 12 This is a specific embodiment of the invention showing the CMP gather and root mean square velocity field of the Inline278 line offset;

[0067] Figure 13 This is a comparison chart showing the effects of single-line and volume offset imaging of Inline278 lines and commercial processing software in a specific embodiment of the present invention. Detailed Implementation

[0068] It should be noted that the following detailed descriptions are exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0069] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments of the present invention. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, and / or combinations thereof.

[0070] like Figure 1 As shown, Figure 1A flow chart of the three-dimensional Kirchhoff integral prestack time migration fast imaging method of the present application, which comprises the following steps:

[0071] Step 100, input the velocity model and seismic data for migration imaging. The seismic data needs to be processed through a series of processes, including noise suppression, multiple wave attenuation, amplitude energy compensation, etc., to prepare a good input trace set for the next step of fast imaging, which can be a common shot trace set or a common midpoint trace set; in addition, the velocity model required in the imaging process is constructed, which can be obtained by using the stacking velocity as the initial migration velocity after velocity analysis, and then obtaining the migration velocity model for time migration after migration velocity analysis; the seismic trace set is evenly distributed to different CPU threads, and the data preparation work for fast time migration is completed.

[0072] Step 101, the travel time of different layers of scattering points on different threads is calculated. In the Kirchhoff integral prestack fast imaging process, the underground formation is regarded as composed of many scattering points. For each scattering point, the travel time is calculated by using the double square root equation, and the wave is assumed to propagate along a straight line from the source point to the scattering point, and the ray path is shown in Figure 2 , S is the source point, R is the receiving point, the total travel time t is the travel time t s from the source point to the scattering point and the travel time t r from the scattering point to the receiving point, and the formula is as follows:

[0073] t=t s +t r

[0074] First, assume that the velocity V is constant, and use the geometric relationship shown in Figure 2 to expand the above formula to obtain the double square root equation:

[0075]

[0076] Where z0 is the depth of the scattering point, x is the position of the midpoint (MP) of the offset relative to the scattering point (SP) located at x=0, and h is one half of the offset distance from the source point to the receiving point. In order to approximately handle the case where the velocity changes strongly with depth and weakly with lateral direction, the double square root equation can be changed to:

[0077]

[0078] Where the migration velocity v mig is the approximate root mean square velocity estimated by Taner and Koehler at t0, and the time t0=t(x=0, h=0) is the average velocity v aveCalculated two-way zero offset time:

[0079]

[0080] In theory, time migration needs to calculate the travel time of each underground imaging grid point, and then pick up the amplitude in the record for stacking imaging. The process of travel time calculation is very complex. If the travel time of each grid point is calculated, a large amount of calculation is needed. However, since the migration velocity field is often a smooth root mean square velocity field, in order to improve the calculation efficiency of the algorithm during the calculation process, not every longitudinal sampling point is calculated by the double square root equation, but according to a certain sampling interval, the coarse imaging grid points are first divided, and then the seismic travel time is calculated on the coarse grid by using the double square root formula, and then the travel time of the intermediate fine grid points is calculated by using linear interpolation or cubic convolution interpolation. The calculation process is shown in the diagram as Figure 3 The basis for this is that the time domain velocity often changes gently, so that by combining coarse and fine grid to calculate travel time, the sudden change of travel time field will not cause calculation error, resulting in inaccurate imaging. By combining coarse and fine grid to calculate travel time, the calculation efficiency can be greatly improved.

[0081] The selection of longitudinal sampling interval△t is an important parameter in fast imaging. Selecting a larger△t can improve the migration calculation efficiency, but the imaging accuracy will be lost. Selecting a smaller△t can ensure the imaging accuracy, but it will also reduce the calculation efficiency. The selection of△t is related to the main frequency f0 of seismic wave. Under the premise that the longitudinal velocity does not change dramatically, the following formula can be used to select△t interval, which can basically ensure the good balance of calculation efficiency and imaging accuracy:

[0082] 2 / f0≤△t≤4 / f0

[0083] In the interval of two△t, the travel time of the sampling point is calculated by using the double square root equation, and in the interval, the travel time is calculated by using the idea of sampling point interpolation. By combining the two, the calculation efficiency of travel time is improved.

[0084] Step 102, in order to further improve the calculation efficiency of the fast imaging algorithm, Streaming SIMD Extensions (SSE) technology is introduced to improve the acceleration ability of the single workstation only relying on CPU calculation on the spot, and the multi-thread parallel technology is used to further improve the algorithm calculation efficiency and meet the needs of the spot fast imaging. Multi-thread calculation can fully exert the calculation advantages of multi-core CPU. Since the threads communicate in the local machine and share the storage mode, the efficiency is very high. The principle of this technology is that the processor with Intel SSE instruction set support has 8 128-bit registers, and each register can store 4 (32-bit) single-precision floating-point numbers. SSE provides an instruction set, and the instructions in the instruction set can allow floating-point numbers to be loaded into these 128-bit registers, and these numbers can be subjected to arithmetic logic operations in these registers. The problem of multi-thread parallel is that the shared storage space cannot open a larger memory space for each thread, and the pre-stack time migration imaging technology developed by the application can well avoid this problem. Since Kirchhoff migration is a one-to-many data mapping process, each imaging line in the three-dimensional imaging space can be sequentially migrated according to a single channel. At this time, the memory overhead does not affect the calculation efficiency. The different threads share the input seismic data, which can effectively improve the calculation performance of multi-core CPU in the pre-stack time migration process. After using the SSE technology, the algorithm based on double square root travel time calculation can be written as follows.

[0085] For every 4 elements in the array (the 4 elements here refer to the travel time of the scattering points at different sampling depths in the vertical direction)

[0086] {

[0087] Load the 4 vertical sampling points in the array into a 128-bit SSE register

[0088] Complete the operation of calculating the double square root travel time of the 4 sampling numbers in one CPU instruction execution period

[0089] Take out the obtained 4 results and write them into the memory

[0090] }

[0091] Without considering the overhead of loading floating-point data in the SSE register, the algorithm can calculate the travel time of 4 vertical sampling points at a time in one CPU period, so the calculation efficiency can be improved by about 4 times.

[0092] Step 103: For each scattering point, perform geometric diffusion compensation for amplitude energy, complete the scanning and superposition of diffraction energy, and realize migration imaging processing based on scattering points. During wave propagation, if scattering and absorption are not considered, when a wave passes through a medium, its initial energy expands on the expanding wavefront, causing a decrease in wave amplitude. This phenomenon is called geometric diffusion. Therefore, eliminating geometric diffusion is a key step in realizing time migration processing based on the integral method. Generally, geometric diffusion compensation to obtain relatively amplitude-preserving imaging results is achieved through the construction and application of amplitude weighting functions. This invention improves and optimizes the weighting function based on the three-dimensional amplitude-preserving weighting function of Bleistein et al. (2001). Based on the derived weighting function formulas for two-dimensional and three-dimensional Kirchhoff pre-stack time migration, compensation for energy loss due to geometric diffusion at scattering points is achieved, improving the amplitude preservation performance of the fast imaging algorithm.

[0093] For two dimensions,

[0094]

[0095] For three dimensions,

[0096]

[0097] Q m Let z represent the weighting function, z represent the depth at the imaging point, and v represent the depth at the imaging point. mig t represents the offset imaging velocity. shot When representing the travel time from the earthquake source to the scattering point, t rec This indicates the travel time from the detector point to the scattering point.

[0098] Step 104, in the three-dimensional case, based on Figure 4 The offset aperture shown completes the scattering at each scattering point ψ. xy The diffracted hyperboloid energy superposition is used to achieve the offset and repositioning of each scattering point. The offset aperture is a cone with its vertex above the ground surface, defined as follows: ψ xy The center position of the scattering point is indicated by θ, and the apex angle of the cone is indicated by θ. Aperture calculation follows:

[0099]

[0100] To achieve this, where v mig t0 represents the offset speed, and t0 represents the round trip time.

[0101] The imaging results from n threads are merged to complete the fast imaging processing of the entire data volume.

[0102] The following are several specific embodiments of the application of this invention:

[0103] Example 1

[0104] In the embodiment 1 of the present application, the Kirchhoff pre-stack time migration fast imaging provided by the present application comprises the following steps:

[0105] Step 100, reading in CMP gathers of 3D seismic data. Taking the actual data of the eastern victory area as an example, the test data is CMP gather data, the CMP number is from 310 to 2165, the CMP interval is 12.5m, the longitudinal time sampling point is 1251, the sampling rate is 4ms, there are 301-inline310 ten lines involved in operation, the data body size is 10.3G, and the line interval between lines is 25m. Taking the inline305 line as an example, the depth domain layer velocity field is as shown in Figure 5 (a), and the extracted CMP gather seismic record is as shown in Figure 5 (b). The data is distributed to 8 threads.

[0106] Step 101, for each data, the travel time at each grid point is calculated every 20-32 sampling points based on the following double square root equation, and the intermediate fine grid point sampling points are calculated by using linear interpolation or cubic convolution interpolation.

[0107]

[0108] Where the migration velocity v mig is the approximate root mean square velocity estimated by Taner and Koehler at t0, and the time t0=t(x=0, h=0) is the two-way zero offset time calculated by the average velocity v ave .

[0109]

[0110] Step 102, in order to further improve the calculation efficiency of the fast imaging algorithm, the single instruction multiple data stream expansion (SSE, Streaming SIMD Extensions) technology is introduced, the SSE instruction set is used in parallel, 4 sample points are loaded into a 128-bit register at a time, the operation of calculating the double square root travel time of the 4 sampling points is completed in a CPU instruction execution period, and then the obtained 4 results are taken out and written into the memory. Through this technology, the calculation efficiency can be improved by about 4 times without considering the overhead of loading floating point data in the SSE register.

[0111] Step 103, for each scattering point, the geometric diffusion compensation of amplitude energy is carried out, the scanning superposition of diffraction energy is completed, and the scattering point based migration imaging processing is realized. Based on the weighted function formula of three-dimensional Kirchhoff pre-stack time migration, the compensation of geometric diffusion loss energy of the scattering point is realized, and the amplitude preservation of the fast imaging algorithm is improved.

[0112] For three-dimensional,

[0113]

[0114] Where Q m represents the weight function, z represents the depth at the imaging point, v mig represents the migration imaging velocity, t shot represents the travel time of the source to the scattering point, t rec represents the travel time of the receiver to the scattering point.

[0115] Step 104, in the three-dimensional case, based on the conical migration aperture, the diffraction hyperboloid energy superposition of each scattering point ψ xy is completed, and the migration homing processing of each scattering point is realized. The migration aperture is a circular cone with a vertex above the ground surface, and its definition is as follows: ψ xy represents the center position of the scattering point, and θ represents the top angle of the circular cone. The aperture calculation is as follows:

[0116]

[0117] is realized, where v mig represents the migration velocity, and t0 represents the two-way travel time.

[0118] Figure 6 is a comparison chart of the migration effects of the rectangular aperture and the conical aperture of the present application, and from the comparison chart, it can be seen that the imaging data signal-to-noise ratio of the conical aperture is better than that of the rectangular aperture because the data contribution within the Fresnel zone radius is considered.

[0119] The imaging results of 8 threads are combined to complete the fast imaging processing of the entire data body.

[0120] Figure 7 is the fast migration imaging result of the inline 305 line, Figure 8is the result of inline 305 line commercial software migration imaging, from the comparison of two figures, it can be seen that the imaging effect of the rapid imaging method developed by the application is equivalent to that of the commercial software imaging, but the time of the rapid imaging is 730s (about 12 minutes), while the migration imaging time of the commercial software Kirchhoff integral method is 2930s, and the calculation efficiency is increased by about 4 times. Table 1 compares the calculation efficiency of imaging a single line in different imaging calculation modes, and from the comparison, it can be seen that after using single machine 8-core parallel plus instruction set parallel, the calculation efficiency is increased by about 14 times compared with single machine single-core serial calculation, so that the rapid imaging can be applied to field site quality monitoring, and the purpose of rapid evaluation of field site data quality is better achieved, and the quality monitoring precision of field site is improved.

[0121] After the single line completes the rapid migration imaging test, ten lines are implemented as a whole to realize the volume migration imaging, the calculation precision and calculation efficiency of the rapid imaging method processing three-dimensional data volume migration are tested. From the imaging effect, it can be seen that since the data contribution in the migration aperture is considered, the data imaging result is better, and the signal-to-noise ratio is higher. Still taking the inline 305 line as an example, the single line migration, adjacent 3-line volume migration and 10-line volume migration of the line are successively performed, and the imaging effects are shown in Figure 9 (a)-(c), respectively. From the comparison figure of the line, it can be seen that the more data used in migration, the greater the data volume contribution in the migration aperture range, and the better the imaging effect. From the calculation efficiency, the data whole rapid imaging takes 142 minutes, while the imaging time of the commercial software under the same configuration is 201 minutes, and the calculation efficiency of the rapid imaging is increased by about 30%. The high efficiency of the calculation efficiency ensures the effectiveness and adaptability of the rapid imaging technology in field site data quality control.

[0122] Embodiment 2

[0123] In the specific embodiment 2 of the application, the trial data is CMP gather data of Luo Jia Inline 350 line, the whole data volume is 2.8G, the data longitudinal time sampling point is 1500, and the sampling rate is 4ms. The test server is 12 threads, the memory is 64G, and the hard disk capacity is 10T. Figure 10 is the single shot record ( Figure 10 a) and the velocity field ( Figure 10 b) of the data, Figure 11 is the result of rapid imaging ( Figure 11 a) and the commercial software ( Figure 11b) The comparison of the two results shows that the imaging effects of both are comparable for depicting medium-deep strata and faults. The rapid imaging effect in local areas is slightly worse. However, the rapid imaging time of this invention is 2460s (about 40 minutes), while the imaging time of commercial software is 60 minutes, which improves the calculation efficiency by 30%.

[0124] Example 3:

[0125] In a specific embodiment 3 of the present invention, the test data is three-dimensional block data of the Binzhou area of ​​Shengli Oilfield, with a total data volume of 13G and 10 survey lines from inline278 to inline287. The device used for testing rapid imaging is a single microcomputer with 16G of memory and 1T of hard disk. Figure 12 This is the single-shot record of the data. Figure 12 a) and velocity field ( Figure 12 b), Figure 13 This is the result of inline278 single-line fast imaging. Figure 13 a) Rapid imaging results of 10 survey lines ( Figure 13 b) Imaging of 10 survey lines using commercial software ( Figure 13 c) Trial calculation results, from Figure 13 The comparison results show that the fast imaging results of the data volume ( Figure 13 b) is better than single-line imaging results ( Figure 13 a) It achieves higher imaging accuracy and better imaging quality for medium-deep strata and faults, and also has a higher signal-to-noise ratio for the entire profile. Compared with the volumetric imaging of commercial software, the imaging effect of both is comparable for medium-deep strata and faults, but the signal-to-noise ratio is slightly worse. The entire migration process takes 4 minutes for single-line fast imaging migration, 42 minutes for fast volumetric imaging migration of this invention, and 120 minutes for volumetric migration imaging of commercial software, which improves the computational efficiency by 35%.

[0126] Finally, it should be noted that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0127] Except for the technical features described in the specification, all other technologies are known to those skilled in the art.

Claims

1. A three-dimensional Kirchhoff integral method for pre-stack time-migrating fast imaging, characterized in that, This three-dimensional Kirchhoff integral pre-stack time-migrating fast imaging method includes: Step 1: Acquire velocity models and seismic data for migration imaging; Step 2: Complete the calculation of the travel time of scattering points at different layers on different threads; Step 3: Calculate the trip of 4 longitudinal sampling points per track in parallel based on the instruction set and store the result in the SSE register; Step 4: For each scattering point, perform geometric diffusion compensation of amplitude energy to complete the scanning and superposition of diffraction energy, and realize offset imaging processing based on scattering points; Step 5: Perform offset imaging within the offset aperture range to complete multi-threaded fast imaging processing in single-machine mode; In step 4, the weighted function formula of the three-dimensional Kirchhoff pre-stack time migration is used to compensate for the energy loss due to geometric diffusion of scattering points, thereby improving the amplitude preservation of the fast imaging algorithm. In step 4, for two dimensions, For three dimensions, Q m Let z represent the weighting function, and z represent the depth of the imaged location at a scattering point, in meters; v mig t represents the offset imaging velocity. shot When representing the travel time from the source point to the scattering point, t rec This indicates the travel time from the detector point to the scattering point.

2. The three-dimensional Kirchhoff integral method for pre-stack time migration fast imaging according to claim 1, characterized in that, In step 1, the seismic data is processed, including noise suppression, multiple attenuation, and amplitude energy compensation, to prepare a good input gather for the next step of rapid imaging. This gather is either a common shot gather or a common center gather.

3. The three-dimensional Kirchhoff integral method for pre-stack time migration fast imaging according to claim 2, characterized in that, In step 1, a velocity model required for the imaging process is constructed. This model uses the stacking velocity obtained after velocity analysis as the initial migration imaging velocity. After migration imaging velocity analysis, a migration imaging velocity model for time migration is obtained. Seismic gathers are then evenly distributed to different CPU threads to complete the data preparation work for fast time migration.

4. The three-dimensional Kirchhoff integral method for pre-stack time migration fast imaging according to claim 1, characterized in that, In step 2, during the Kirchhoff integral method pre-stack rapid imaging process, the subsurface strata are considered to consist of many scattering points. For each scattering point, the travel time is calculated using a double square root equation. During migration, it is assumed that the wave propagates in a straight line from the source point to the scattering point. The total travel time t is the travel time t from the source point to the scattering point. shot and the travel time t from the scattering point to the receiving point rec The sum, the formula is as follows: t=t shot +t rec First, assuming the velocity υ is constant, expanding the above equation yields the double square root equation: Where: z0 is the depth of the scattering point, in meters; x is the position of the midpoint of the source-receiver distance MP relative to the scattering point SP located at x=0; h is half the offset distance from the source point to the receiver point; and υ is the velocity. To approximate the case where the velocity varies strongly with depth but weakly with lateral movement, the double square root equation is modified as follows: Where t0 is the average velocity v ave The calculated two-way zero-shot-receiver distance time and the offset imaging velocity v mig The approximate root mean square velocity estimated by Taner and Koehler at t0 is:

5. The three-dimensional Kirchhoff integral method for pre-stack time migration fast imaging according to claim 4, characterized in that, In step 2, when calculating the time offset, according to a certain sampling interval Δt, the coarse imaging grid points are first divided, and the seismic travel time is calculated on the coarse grid using the double square root formula. Then, the travel time of the intermediate fine grid points is obtained by using linear interpolation or cubic convolution interpolation.

6. The three-dimensional Kirchhoff integral method for pre-stack time migration fast imaging according to claim 5, characterized in that, In step 2, the selection of Δt is related to the dominant frequency f0 of the seismic wave. Under the premise that the longitudinal velocity does not change drastically, the following formula is used to select the Δt interval, which can basically ensure a good balance between computational efficiency and imaging accuracy: 2 / f0≤△t≤4 / f0 At the interval of two Δt, the travel time of the sample points is calculated using the double square root equation, while within the interval, the travel time is calculated using the idea of ​​sample point interpolation. By combining the two methods, the efficiency of travel time calculation is improved.

7. The three-dimensional Kirchhoff integral method for pre-stack time migration fast imaging according to claim 1, characterized in that, In step 3, under the condition that there is only a single workstation on site, the calculation and storage are completed by relying on the acceleration capability of the CPU. The Single Instruction Multiple Data Stream Extended (SSE) technology is introduced. The processor with Intel SSE instruction set support has eight 128-bit registers, each of which can store four 32-bit single-precision floating-point numbers. These numbers can then be used for arithmetic and logical operations in these registers.

8. The three-dimensional Kirchhoff integral method for pre-stack time migration fast imaging according to claim 7, characterized in that, In step 3, after using the SSE technique, the algorithm for calculating travel time based on the double square root is as follows: For each of the four elements in the array, where the four elements refer to the scattering points at different vertical sampling depths during the travel, these four vertical sampling points in the array are loaded into a 128-bit SSE register; the double square root travel operation of these four sample numbers is completed in one CPU instruction execution cycle; and the four results are retrieved and written into memory. Without considering the overhead of loading floating-point data into the SSE register, the algorithm can improve computational efficiency by about 4 times when calculating the travel of 4 vertical sampling points in one CPU cycle; by combining instruction set parallelism, the on-site rapid imaging capability of a single workstation is improved.

Images