A method for calculating a three-dimensional seismic imaging point spread function

Abstract

The present application relates to a kind of three-dimensional seismic imaging point spread function calculation method, belong to the technical field of seismic exploration.The method disperses the huge five-dimensional travel time memory in all nodes, and the travel time slice required by current target line is reduced into the current computing point memory.The present application makes full use of the computing and storage capacity of current multi-node cluster, matches the appropriate parallel framework, and can efficiently calculate three-dimensional point spread function.Using the calculated point spread function for imaging amplitude correction and least square migration can improve the reliability of imaging results.

Description

TECHNICAL FIELD

[0001] The present application relates to a three-dimensional seismic imaging point spread function calculation method, belonging to the technical field of seismic exploration. BACKGROUND

[0002] Migration is one of the core technologies of seismic data processing in oil and gas exploration. The goal of traditional migration techniques is to locate subsurface structures, while inversion imaging represented by least squares migration can better estimate the subsurface true reflection coefficients. The book Seismic Inversion by Schuster in 2017 describes the method principle of least squares migration and the challenges in method theory and implementation. In the 2019 SEG annual conference paper Leastsquares migration: Current and future directions, Schuster pointed out that the computational cost of least squares migration in data domain and imaging domain is very large, especially in data domain, which is 10-30 times the computational cost, and there is also the problem of convergence. In the 2005 Geophysics article Directional illumination analysis using beamlet decomposition and propagation by Ru-Shan Wu et al. and the 2006 Geophysics article Wave-equation based seismic illumination analysis by Xiao-Bi Xie et al., it is pointed out that from the perspective of imaging correction of illumination, the subsurface true reflection coefficient can also be obtained, and it is more stable, and the difficulty lies in the resolution function representing the subsurface illumination and propagation factors, that is, the computational cost of point spread function is increased. Subsequently, Chen and Xie proposed an efficient algorithm for calculating two-dimensional point spread function in the 2015 SEG annual conference paper An efficient method for broadband seismic illumination and resolution analyses. In the 2019 paper of Oil and Gas Exploration based on point spread function inversion imaging of viscous medium by Wang Xuejun and Ren Haoran et al., the point spread function is calculated to carry out two-dimensional viscous medium inversion imaging.

[0003] Currently, two-dimensional and three-dimensional migration imaging techniques are relatively mature and can be applied to industrial production. The calculation amount of seismic inversion imaging represented by three-dimensional least square migration is still challenging, and the calculation and storage of a full Hessian matrix in three-dimensional seismic inversion are almost impossible (for example, a small three-dimensional model with a volume of 1G, the size of the full Hessian is the square of the model, that is, 1T, and the memory of the current computer cannot meet the requirement). SUMMARY

[0004] The purpose of the present application is to provide a three-dimensional seismic imaging point spread function calculation method to solve the problems of large calculation amount and low efficiency of the current calculation method.

[0005] The present application provides a three-dimensional seismic imaging point spread function calculation method to solve the above technical problems, which comprises the following steps:

[0006] 1) Pile point arrangement is performed within the range of the observation system, and each pile point is arranged on the surface of the observation system;

[0007] 2) A plurality of target lines are set, and the seismic data volume in the observation system is profiled along any target line;

[0008] 3) The shot-receiver pair information of any target line is obtained, and the number of corresponding travel time field pile points is determined according to the shot-receiver pair range;

[0009] 4) The three-dimensional travel time field of each pile point under the target line is obtained, and parallel calculation is performed through a plurality of calculation nodes to obtain the travel time field of each pile point;

[0010] 5) The PSF data information of the target line is obtained by using the travel time field data of each pile point under the target line; and the PSF of all target lines is integrated to obtain the PSF of the entire observation system.

[0011] The present application sets a plurality of target lines along the set direction of the observation system, and profiles the observation system data volume along each target line; then, the travel time slices required by the current target line are reduced into the current calculation point memory. The present application maximally utilizes the calculation and storage capacity of the current multi-node cluster, matches a proper parallel framework, and can efficiently perform three-dimensional point spread function calculation.

[0012] Further, the target line is set along the X-axis or Y-axis direction of the observation system.

[0013] Further, the step 4) averages the travel time field of the pile points and distributes it to each calculation node according to the number of pile points.

[0014] Further, the memory consumption of a single computing node PSF is PSF (Npsf_str, Npsf / Node), and the memory consumption of a single computing node in the traveltime is Traveltime_slices (Ntz, Ntx, 2, Npile), Node is the number of computing nodes, Npile is the number of pile points, (Ntz, Ntx, 2) is the calculation scale of a single pile point traveltime on the target line, Npsf_str is the size of a point spread function calculation psf structure.

[0015] Further, the traveltime field of each image point of the observation system is saved in all computing nodes, and the pile point traveltime field used by the target line is sent to each computing node through broadcasting at each time of calling the target line.

[0016] Further, the method uses OpenMP for parallel processing.

[0017] Further, each pile point in the step 1) is uniformly arranged on the surface of the observation system. BRIEF DESCRIPTION OF DRAWINGS

[0018] Figure 1 is a flowchart of the three-dimensional seismic imaging point spread function calculation method of the present application;

[0019] Figure 2 is a schematic diagram of the PSF basic concept in the present application;

[0020] Figure 3 is a schematic diagram of the three-dimensional PSF calculation implementation in the present application;

[0021] Figure 4 is a schematic diagram of the PSF of the three-dimensional Salt model Line 300 in the embodiment of the present application;

[0022] Figure 5-a is a schematic diagram of the conventional Kirchhoff integral migration result of the three-dimensional Salt model Line 300 in the embodiment of the present application;

[0023] Figure 5-b is a schematic diagram of the PSF corrected migration result of the three-dimensional Salt model Line 300 in the embodiment of the present application. DETAILED DESCRIPTION

[0024] The specific embodiments of the present application will be further described below in combination with the drawings.

[0025] In classical seismic data acquisition and imaging systems, the Hessian matrix represents the system response encompassing both acquisition (forward propagation of seismic waves) and imaging (backward propagation of seismic waves). The point spread function, on the other hand, is the Hessian matrix of local support and can be viewed as the impulse response of the seismic data acquisition and imaging system. The core issue in seismic data inversion imaging is to correct for systemic effects caused by observation and propagation, ensuring accurate imaging of subsurface parameters. In reflection seismic exploration, the seismic data acquisition or forward modeling process is as follows:

[0026] Lm = d obs (1)

[0027] Where L represents the wave propagation operator, m represents the subsurface reflection parameter, and d obs This represents the collected data. The exact solution for the subsurface reflection parameter *m* can be obtained by solving the inverse of the propagation operator. Typically, the inverse of *L* does not exist. To estimate the subsurface reflection parameter, migration is introduced in exploration seismic surveys. The observed data is directly migrated back-propagated to obtain the structural information of the subsurface reflection parameter. The process is as follows:

[0028]

[0029] in, Structural imaging results of subsurface reflection parameters, L T This represents the migration operator (the transpose of the forward propagation operator). However, the classical migration process described above only yields a relatively blurry reflection structure of the subsurface parameters. More accurate estimates of the subsurface parameters can be obtained through iterative inversion least-squares migration. The target functional for inversion imaging in the least-squares sense can be written as...

[0030]

[0031] Generally, the above objective function can be solved using gradient-guided methods. The normal equation for the above formula is as follows:

[0032]

[0033] Where H = L T L is the second derivative of the objective function in equation (3), which is the Hessian matrix of the linear operator L. From the above equation, it can be seen that the result of conventional imaging is equivalent to a blurred result of the Hessian matrix applied to the subsurface true reflectance coefficient. The imaging domain least-squares migration attempts to directly solve the normal equation (4), estimate the inverse of the Hessian matrix, and apply it to the conventional imaging result, thereby obtaining the inverted imaging result.

[0034]

[0035] in, Hessian corrected comparison There are higher precision least squares migration results.

[0036] In the Born approximation, the positive operator in equation (1) is explicitly written as

[0037] d obs (x r , x s , ω) = ∫ω 2 G(x, x r , ω)G(x, x s , ω) s(ω) m(x) dV (6)

[0038] In the above equation, d obs (x r , x s , ω) represents the seismic acquisition (forward) data with the surface shotpoint at (x r , x s ) and the frequency domain, s(ω) represents the wavelet, and G(x, x s , ω) and G(x, x r , ω) are the Green's functions of the shotpoint and the receiver, respectively. By substituting the above equation into the Hessian matrix, H = L T L, we can obtain

[0039]

[0040] where ξ = (x s , x r ) represents the surface shot-receiver pair, and x i and x j both represent a grid point in the subsurface. Based on the asymptotic approximation theory, the Green's function under the WKBJ approximation can be expressed as

[0041]

[0042] where T(x, x s ) represents the travel time term from the shot (receiver) point to the subsurface imaging point, and A(x, x s ) represents the corresponding amplitude term. By substituting equation (8) into equation (7), the expression of the Hessian matrix can be obtained as follows

[0043]

[0044] where R ss is the cross-correlation function of the wavelet. In the above equation, the dimensions of x i and x j are both the model space, and it can be seen that the size of the Hessian matrix is the square of the model. However, the Hessian matrix is sparse and strip-shaped, which means that its elements have the characteristics of local support, that is, they are related to only a single point x iOnly its surrounding x j is valid. According to its locality, for any point x i , only a limited range of x j around it is calculated as the point spread function of the point.

[0045] According to formula (9), in the three-dimensional case, the calculation amount and storage amount of calculating all the point spread functions of the local support are very large, and the difficulty in the implementation of the three-dimensional point spread function calculation is the management of the five-dimensional travel time and amplitude field. Therefore, the application provides a three-dimensional seismic imaging point spread function calculation method, which adopts a target line output-based mode, only extracts and the travel time and amplitude related to the target line are included in the current calculation node, and the vast amount of travel time and amplitude unrelated to the target line is stored in the memory of all calculation nodes; each calculation node adopts a parallel calculation mode, and the parallel structure thereof is as follows: the outermost loop of the parallel framework according to the target line output is the target line loop, the travel time and amplitude field related to the target line is shared in the loop, and the data projection in the observation system related to the target line is processed. The implementation process of the method is shown in Figure 1 , and the specific implementation process is as follows.

[0046] 1. Obtain an observation system file, and arrange pile points in the observation system range.

[0047] The application arranges a plurality of pile points in the observation system range, performs parallel travel time calculation by using Nnode calculation nodes, and each calculation node can be allocated Npile / Nnode pile point travel time calculation. Npile pile points are evenly arranged on the surface of the observation system, and each pile point corresponds to a certain number of image points. The size of single pile point travel time calculation is (Ntz, Ntx, Nty), and the Npile / Nnode pile point travel time allocated to the calculation node can be parallel by using multiple processes. After the travel time calculation is parallel, the travel time field obtained by the single calculation node has a size of (Ntz, Ntx, Nty, Npile / Nnode), and all calculation nodes can obtain the full five-dimensional travel time field (Npile is the two-dimensional pile point on the ground).

[0048] 2. Set a plurality of target lines along the X-axis or Y-axis direction of the observation system, and divide the observation system along the target lines.

[0049] This invention allows for target line setting based on actual conditions. The target line can be set along the X-axis of the observation system or along the Y-axis. In this embodiment, the output target line coordinates are used as a reference. Two Y-axis slices of travel time fields near the target line are extracted. Based on the extracted travel time fields, the travel time fields on all computing nodes are reduced. The reduced travel time field for each computing node is Traveltime_slices(Ntz,Ntx,2,Npile / Nnode).

[0050] 3. Perform parallel processing on the NPSF within the target line.

[0051] After dividing the target line, each profile on both sides of the target line contains several stake points. Parallel processing of the NPSF within the target line is required. The spread function calculation for a single stake point is called the PSF, such as... Figure 2 As shown, the boxed portion represents the calculation range for the travel time and amplitude of the corresponding pile point, while the center point of the PSF and its neighborhood coordinates are x and x respectively. i and x j A point-like diffusion function computation is called a PSF. The size of each PSF structure is nx_psf * ny_psf * nz_psf = Npsf_str. The memory overhead of a single compute node's PSF is PSF(Npsf_str, Npsf / Node). The travel time memory required by a single compute node is Traveltime_slices(Ntz, Ntx, 2, Npile). The PSF structure and PSF index file are written in parallel.

[0052] A schematic diagram of 3D PSF calculation is shown below. Figure 3 As shown, the travel time field (Ntz, Ntx, Nty, Npile) represented by the cube on the right is always stored in the memory of N compute nodes. Each time it is called, the travel time Traveltime_slices(Ntz, Ntx, 2, Npile) shown in the profile slice is reduced and broadcast to all compute nodes.

[0053] The process of PSF calculation is as follows: Figure 1 As shown, it is assumed that the observation system targeted by this invention has N lines along a set direction. Starting from target line 1, the shot-receiver pair information Nsr of the target line is first obtained, and the number of corresponding travel time field stakes is determined according to the range of the shot-receiver pair; the underground three-dimensional travel time field of each stake is obtained; the total travel time field of the target line is synchronously extracted by calculating the node reduction, and the PSF data information of the target line is obtained; following the above process, the remaining target lines are also processed to obtain the PSF of each target line; by combining the PSF of all target lines, the PSF of the entire observation system can be obtained.

[0054] Through the above process, the calculation of the three-dimensional seismic imaging point spread function can be realized. Taking the three-dimensional Salt model Line 300 in the embodiment as an example, the calculation result of the PSF is as shown in the figure Figure 4 The imaging amplitude correction and the least square migration are performed by using the calculation result. The conventional Kirchhoff integral migration result of the three-dimensional Salt model Line 300 is as shown in the figure Figure 5-a The PSF corrected migration result is as shown in the figure Figure 5-b It can be known by comparison that the reliability of the image obtained by the spread function calculation and correction according to the application is higher.

[0055] The application is aimed at the current characteristics of multiple computer nodes, independent memory in a single node and large total memory. The huge five-dimensional travel time memory is dispersedly placed in all the computing nodes, and the travel time slice required by the current target line is reduced into the current computing node memory. The application maximally utilizes the computing and storage capacity of the current multi-node cluster, matches a proper parallel framework, and can efficiently perform the calculation of the three-dimensional point spread function. The imaging amplitude correction and the least square migration are performed by using the calculated point spread function, and the reliability of the imaging result can be improved.

Claims

1. A method for three-dimensional seismic imaging point spread function computation, characterized by, The method comprises the following steps: 1) arranging pile points in the range of the observation system, each pile point being arranged on the surface of the observation system; 2) arranging a plurality of target lines, and profiling the seismic data volume in the observation system along any target line; 3) obtaining the shot-receiver pair information of any target line, and determining the number of corresponding travel time field pile points according to the shot-receiver pair range; 4) obtaining the three-dimensional underground travel time field of each pile point under the target line, and performing parallel calculation through a plurality of calculation nodes to obtain the travel time field of each pile point; 5) obtaining the PSF data information of the target line by using the travel time field data of each pile point under the target line; and integrating the PSFs of all target lines to obtain the PSF of the entire observation system. The memory consumption of a single calculation node PSF is PSF (Npsf_str, Npsf / Node), and the travel time memory to be called by a single calculation node is Traveltime_slices (Ntz, Ntx, 2, Npile), Node is the number of calculation nodes, Npile is the number of pile points, (Ntz, Ntx, 2) is the calculation scale of a single pile point travel time on a target line, and Npsf_str is the size of a pixel spread function calculation psf structure.

2. The three-dimensional seismic imaging point spread function calculation method according to claim 1, characterized in that: the target line is arranged along the X-axis or Y-axis direction of the observation system.

3. The method of claim 1, wherein, In step 4), the travel time field of the pile point is evenly distributed to each calculation node according to the number of pile points.

4. The method of claim 1, wherein, The travel time field of each pixel of the observation system is stored in all calculation nodes, and the travel time field of the pile point used by the target line is sent to each calculation node through broadcasting when the target line is called each time.

5. The method of claim 1, wherein, OpenMP is used for parallel processing.

6. The method of claim 1, wherein, In step 1), each pile point is uniformly arranged on the surface of the observation system.

Images