Cross-node high-parallelism point spread function acquisition method and system

Abstract

The invention provides a cross-node high-parallel point spread function acquisition method and system, and belongs to the technical field of seismic wave imaging. The method comprises the following steps: determining scattering points according to a reflection coefficient model; then different scattering points are distributed to different nodes; different offset pairs corresponding to each scattering point are distributed to different calculation cores on the same node for distributed parallel calculation; and finally, obtaining a local point spread function through a calculation result of each node, and obtaining a global point spread function according to all the local point spread functions. According to the method, the angle and the amplitude field on the ray path are solved in a ray tracing mode, and the offset pairs corresponding to different scattering points are distributed to different calculation cores of different nodes, so that the calculation efficiency of the point spread function is greatly improved, the subsequent application of technologies such as rapid imaging is facilitated, and the method is suitable for large-scale popularization and application. And further development of technologies such as rapid imaging in an actual work area is promoted.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic wave imaging, and particularly relates to a method and system for obtaining a point spread function with high cross-node parallelism. Background Art

[0002] Mathematically speaking, the point spread function can be represented by a row (or a column) of elements in the Hessian matrix and can be used as an approximation of the Hessian matrix. At present, the point spread function has been widely used in seismic exploration. In the least-squares imaging for obtaining the reflection coefficient through the point spread function and the conventional profile, the depth-domain interpretation of inverting the subsurface impedance through the point spread function, or the fast imaging algorithm through the point spread function and the reflection coefficient profile, the point spread function plays an important bridging role.

[0003] However, in actual production, the currently commonly used scheme for calculating the point spread function by the wave equation far fails to meet the requirements of actual production. Among them, the fast imaging technology mainly analyzes the imaging profiles under different observation systems to determine whether the observation system parameters are appropriate. Therefore, frequent calculations of the point spread function are required. The original scheme of calculating by the wave equation has a long calculation time and cannot meet the production requirements. Summary of the Invention

[0004] The purpose of the present invention is to solve the problems existing in the above-mentioned prior art, and provide a method and system for obtaining a point spread function with high cross-node parallelism. By adopting the ray tracing method to obtain the angles and amplitudes on the ray path, and at the same time adopting a high cross-node parallel computing strategy, the fast calculation of the point spread function is realized, and the calculation efficiency is improved.

[0005] The present invention is realized by the following technical solutions:

[0006] In the first aspect of the present invention, a method for obtaining a point spread function with high cross-node parallelism is provided. The method first determines scattering points according to the reflection coefficient model; then distributes different scattering points to different nodes; then distributes different shot-receiver pairs corresponding to each scattering point to different computing cores on the same node for distributed parallel computing; finally, obtains local point spread functions through the calculation results of each node, and obtains the global point spread function according to all local point spread functions.

[0007] A further improvement of the present invention lies in:

[0008] The method includes:

[0009] (1) Input the observation system, background velocity, and seismic wavelet;

[0010] (2) Determine each scattering point in the point spread function model according to the reflection coefficient model;

[0011] (3) Assign different scattering points in the point spread function model to different nodes;

[0012] (4) Obtain the local point spread function corresponding to each scattering point;

[0013] (5) Obtain the global point spread function;

[0014] (6) Output the global point spread function.

[0015] A further improvement of the present invention lies in:

[0016] The operation of determining each scattering point in the point spread function model in step (2) includes:

[0017] First, start searching from the first point in the reflection coefficient model, search for the first point with a non-zero value, and set this point as the first scattering point;

[0018] Then, search for the next point with a non-zero value at a set interval, and set these points as the second scattering point, the third scattering point until the last scattering point in sequence;

[0019] Finally, set the values of all scattering points in the reflection coefficient model to 1, and set the values of other points to 0 to obtain the point spread function model.

[0020] A further improvement of the present invention lies in:

[0021] The operation of step (3) includes:

[0022] Assign different scattering points in the point spread function model to different nodes through MPI.

[0023] A further improvement of the present invention lies in:

[0024] The operation of step (4) includes:

[0025] Perform the following operations on each scattering point:

[0026] (41), Assign the shot-receiver pairs corresponding to the scattering point to different computing cores on the same node, and perform distributed parallel computing to obtain the local plane wave corresponding to each shot-receiver pair;

[0027] (42), Superimpose the local plane waves corresponding to all shot-receiver pairs to obtain the local point spread function corresponding to the scattering point.

[0028] A further improvement of the present invention lies in:

[0029] The operation of performing distributed parallel computing in step (41) to obtain the local plane wave corresponding to each shot-receiver pair includes:

[0030] In each computing core, ray tracing is performed starting from the shot point to obtain the travel time field T of the shot point s and the amplitude field A s , and ray tracing is performed starting from the receiver point to obtain the travel time field T of the receiver point r and the amplitude field A r ;

[0031] The angle is calculated using the following formula:

[0032]

[0033]

[0034]

[0035] where P S , P R are the incident directions of the shot point and the receiver point respectively, θ is the angle, and x and z are the directions of the spatial X and Z axes;

[0036] The local plane wave is calculated using the following formula:

[0037] f′ = A s A r ∫||F(w)|| 2 exp(-iw(T s +T r ))+θ)dw

[0038] where f′ is the local plane wave, w is the frequency, and F(w) is the result after Fourier transform of the input seismic wavelet.

[0039] A further improvement of the present invention lies in:

[0040] The operation of step (5) includes: superimposing the local point spread functions corresponding to all scattering points to obtain the global point spread function.

[0041] The second aspect of the present invention provides a system for obtaining a point spread function with high cross-node parallelism, and the system includes:

[0042] An input unit for inputting an observation system, a background velocity, and a seismic wavelet;

[0043] A scattering point determination unit, connected to the input unit, for determining each scattering point in the point spread function model according to the reflection coefficient model;

[0044] An allocation unit, connected to the scattering point determination unit, for allocating different scattering points in the point spread function model to different nodes;

[0045] A local point spread function acquisition unit, connected to the allocation unit, for obtaining the local point spread function corresponding to each scattering point;

[0046] A global point spread function acquisition unit, connected to the local point spread function acquisition unit, for obtaining the global point spread function;

[0047] An output unit, connected to the global point spread function acquisition unit, for outputting the global point spread function.

[0048] In a third aspect of the present invention, there is provided a computer-readable storage medium storing at least one computer-executable program, and when the at least one program is executed by the computer, the computer is caused to execute the steps in the above-mentioned cross-node high-parallel point spread function acquisition method.

[0049] In a fourth aspect of the present invention, there is provided a computer device including a memory and a processor, the memory storing a computer program, and when the computer program is executed by the processor, the processor is caused to execute the steps of the cross-node high-parallel point spread function acquisition method as described above.

[0050] Compared with the prior art, the beneficial effects of the present invention are as follows: The present invention uses ray tracing to obtain the angles and amplitude fields on the ray paths, and allocates the shot-receiver pairs corresponding to different scattering points to different computing cores on different nodes, greatly improving the calculation efficiency of the point spread function, contributing to the subsequent applications of technologies such as fast imaging, and promoting the further development of technologies such as fast imaging in actual work areas. BRIEF DESCRIPTION OF THE DRAWINGS

[0051] Figure 1 is a block diagram of the steps of the method of the present invention;

[0052] Figure 2 is the velocity model in the embodiment;

[0053] Figure 3 is the reflection coefficient model in the embodiment;

[0054] Figure 4 is the observation system diagram in the embodiment;

[0055] Figure 5 is the result of the point spread function for a single node;

[0056] Figure 6 is the result of the point spread function for all;

[0057] Figure 7Results of fast imaging. Detailed implementation mode

[0058] The present invention will be further described in detail below with reference to the accompanying drawings:

[0059] As an important bridge connecting seismic profiles and subsurface reflection coefficients, the point spread function plays an important role in least-squares migration, fast imaging, and depth-domain interpretation. However, in 3D actual data, the calculation speed of algorithms based on wave methods far fails to meet the requirements of timeliness. At the same time, the calculation method of the point spread function for a single node has poor parallelism and low calculation efficiency, which also affects the application of the point spread function in actual work areas. All these seriously hinder the further development of technologies such as fast imaging.

[0060] Aiming at the problems existing in the above-mentioned point spread function, the present invention provides a method for obtaining a cross-node highly parallel point spread function, which mainly realizes the calculation of the point spread function through ray tracing, and further improves the calculation efficiency through a cross-node highly parallel strategy.

[0061] The method of the present invention first determines scattering points according to the reflection coefficient model; then distributes different scattering points to different nodes to achieve parallelism at different scattering point levels; then distributes different shot-receiver pairs corresponding to each scattering point to different calculation cores on the same node for distributed parallel calculation to achieve parallelism at the shot-receiver pair level; finally, obtains the imaging of a single scattering point, that is, the local point spread function, based on the calculation results of each node, and obtains the global point spread function according to all local point spread functions.

[0062] As Figure 1 shown, the method of the present invention includes:

[0063] (1) Input the acquisition system, background velocity, and seismic wavelet; the information of the acquisition system includes: shot point position and geophone position information;

[0064] (2) Reflection coefficient addressing: Determine each scattering point in the point spread function model according to the reflection coefficient model;

[0065] (3) Distribute different scattering points in the point spread function model to different nodes;

[0066] (4) Obtain the local point spread function corresponding to each scattering point;

[0067] (5) Obtain the global point spread function;

[0068] (6) Output the global point spread function.

[0069] Embodiments of the present invention are as follows:

[0070] Embodiment 1:

[0071] The reflection coefficient model in step (2) is obtained as follows:

[0072] Either directly input the reflection coefficient model or, according to the background velocity, obtain the reflection coefficient model according to the formula where v is the velocity, δv is the velocity change, which can be obtained by taking the difference between the input velocity and the result after velocity smoothing.

[0073] The operations for determining each scattering point in the point spread function model in step (2) include:

[0074] First, start searching from the first point in the reflection coefficient model, search for the first point with a non-zero value, set this point as the first scattering point, and the position information of this point is the position information of the first scattering point;

[0075] Then, search for the next point with a non-zero value at a set interval (for example, after finding the first point with a non-zero value, start from this point and search for a point with a non-zero value again after 10 points, and then search for the next point with a non-zero value from this point after another 10 points, and so on until the last point with a non-zero value is found). Set these points as the second scattering point, the third scattering point until the last scattering point in turn. In this way, all scattering points are obtained, and the position information of all scattering points is also determined.

[0076] Finally, set the values of all scattering points in the reflection coefficient model to 1 and the values of other points to 0 to obtain the point spread function model.

[0077] Embodiment 2:

[0078] The operations in step (3) include:

[0079] Allocate different scattering points in the point spread function model to different nodes through MPI (MultiPoint Interface). The nodes and the scattering points do not need to correspond one by one. This is a parallel operation. At the same time, different scattering points are allocated to different nodes. After the calculation of this scattering point is completed, the next scattering point can be allocated to this node.

[0080] Embodiment 3:

[0081] The operations in step (4) include:

[0082] Perform the following operations on each scattering point:

[0083] (41) Allocate the shot-receiver pairs corresponding to the scattering points to different computing cores on the same node, and perform distributed parallel computing to obtain the local plane waves corresponding to each shot-receiver pair.

[0084] Specifically, the operation of "assigning the shot-receiver pairs corresponding to the scattering points to different computing cores of the same node" is automatically performed using existing methods. For example, the built-in multiprocessing module in Python can be used to put the shot-receiver pairs into the process pool of the multiprocessing module, and the multiprocessing module will automatically assign the shot-receiver pairs to different computing cores.

[0085] The operation of "performing distributed parallel computing to obtain the local plane wave corresponding to each shot-receiver pair" includes:

[0086] In each computing core, calculate the angle and amplitude field corresponding to each shot-receiver pair, and then obtain the local plane wave corresponding to each shot-receiver pair through seismic wavelet loading and amplitude correction. The specific calculation method is as follows:

[0087] The travel-time field can be achieved by kinematic ray tracing, and the amplitude field can be achieved by dynamic ray tracing. These are all existing technologies, that is, ray tracing starting from the shot point can obtain the travel-time field T s and amplitude field A s , and ray tracing starting from the receiver point can obtain the travel-time field T r and amplitude field A r , and the calculation of the angle is shown in formula (1):;

[0088]

[0089] where P S , P R are the incident directions of the shot point and the receiver point respectively, θ is the angle (i.e., dip angle), x and z are the spatial X and Z axis directions, represents the solution of the spatial derivative in the X direction, and the same applies to the Z axis.

[0090] Calculate the local plane wave according to formula (2):

[0091] f′ = A s A r ∫||F(w)|| 2 exp(-iw(T s +T r ))+θ)dw (2)

[0092] where f′ is the local plane wave, w is the frequency, and F(w) is the result after Fourier transform of the input seismic wavelet. T s and T r are the travel-time fields of the shot point and the receiver point respectively, θ is the angle, A s and A r are the amplitude fields of the shot point and the receiver point respectively.

[0093] (42), Superimpose the local plane waves corresponding to all shot-receiver pairs to obtain the local point spread function corresponding to the scattering point, that is, superimpose the local plane waves obtained on all calculation cores in step (41) to obtain the local point spread function corresponding to one scattering point.

[0094] In parallel computing, the operation of superimposing the calculation results on different calculation cores (or nodes) is called reduction. In actual processing, the reduction function Reduce can be directly called to implement this superimposing operation. The result obtained after superimposing the plane waves is the point spread function corresponding to the scattering point, which is called the local point spread function.

[0095] The scattering point is the point determined on the reflectivity model, that is, the point whose value is set to 1 in step (1). The shot-receiver pair refers to the combination of a shot point and a geophone point. When a shot point is excited and the seismic wave propagates downward to the scattering point, the geophone point to which the scattering point scatters upward. One shot-receiver pair can be approximated as a plane wave around the scattering point, and the present invention uses this process for calculation.

[0096] Example 4:

[0097] The operation of step (5) includes: superimposing the local point spread functions corresponding to all scattering points to obtain the global point spread function.

[0098] Specifically, the Reduce function can be called again to receive the local point spread functions generated in step (4), that is, reduce the local point spread functions generated by different nodes in step (4) to obtain the final global point spread function.

[0099] Example 5:

[0100] Next, taking a three-dimensional complex velocity model as an example, the present invention will be further described in conjunction with the accompanying drawings and specific implementation manners.

[0101] (1) Read the velocity model v, as Figure 2 shown, simply calculate the reflection coefficient model m ref , as Figure 3 shown;

[0102] (2) Take the information set of the observation system as the input of the function, as Figure 4 shown, and distribute the position information of different scattering points of the point spread function model to different nodes through MPI;

[0103] (3) Automatically distribute the shot-receiver pairs corresponding to each scattering point to different calculation cores on the same node for distributed parallel computing, calculate the angle and amplitude field corresponding to each shot-receiver pair, and form a corresponding plane wave;

[0104] (4) Define a Reduce function that receives a plane wave generated by each shot-receiver pair in step (3) and performs superposition processing to form a local point spread function m corresponding to an underground scattering point locpsf , as Figure 5 shown.

[0105] (5) Then define another Reduce function that receives the local point spread function m generated in step (4) locpsf , and perform reduction processing (using the reduction function MPI_Reduce for reduction, and directly superimposing the results generated by different nodes in step (4) on the master node);

[0106] (6) Output the result of the reduction processing (i.e., the global point spread function) to obtain the final global point spread function m psf , as Figure 6 shown;

[0107] Figure 5 is the point spread function calculated by a single node, Figure 6 is the global point spread function calculated by all nodes. Through the parallel calculation method of scattering points, its calculation efficiency increases linearly with the number of calculation nodes. Figure 7 is the fast imaging result formed by the point spread function of Figure 6 , which indicates that the method proposed by the present invention can well apply the calculated point spread function to fast imaging.

[0108] By distributing scattering points across nodes to form parallel processing at the scattering point level and then performing distributed parallel calculation on the shot-receiver pairs in a single node, the present invention greatly improves the calculation efficiency of the point spread function and promotes the development of subsequent fast imaging technology.

[0109] The above technical solution is only one implementation manner of the present invention. For those skilled in the art, based on the disclosed principle of the present invention, it is very easy to make various types of improvements or deformations, not limited to the technical solution described in the above specific embodiments of the present invention. Therefore, the foregoing description is only preferred and does not have a restrictive meaning.

Claims

1. A method for obtaining a point spread function with high cross - node parallelism, characterized in that: the method first determines scattering points according to a reflection coefficient model; then assigns different scattering points to different nodes; then assigns different shot - receiver pairs corresponding to each scattering point to different computing cores on the same node for distributed parallel computing; finally, obtains a local point spread function from the calculation results of each node, and obtains a global point spread function based on all local point spread functions.

2. The method for obtaining a point spread function with high cross - node parallelism according to claim 1, characterized in that: the method includes: (1) Inputting an observation system, background velocity, and seismic wavelet; (2) Determining each scattering point in the point spread function model according to the reflection coefficient model; (3) Assigning different scattering points in the point spread function model to different nodes; (4) Obtaining a local point spread function corresponding to each scattering point; (5) Obtaining a global point spread function; (6) Outputting the global point spread function.

3. The method for obtaining a point spread function with high cross - node parallelism according to claim 2, characterized in that: the operation of determining each scattering point in the point spread function model in step (2) includes: First, start searching from the first point in the reflection coefficient model, search for the first point with a non - zero value, and set this point as the first scattering point; Then, search for the next point with a non - zero value at a set interval, and set these points as the second scattering point, the third scattering point until the last scattering point in turn; Finally, set the values of all scattering points in the reflection coefficient model to 1, and set the values of other points to 0 to obtain the point spread function model.

4. The method for obtaining a point spread function with high cross - node parallelism according to claim 2, characterized in that: the operation of step (3) includes: Assigning different scattering points in the point spread function model to different nodes through MPI.

5. The method for obtaining a point spread function with high cross - node parallelism according to claim 2, characterized in that: the operation of step (4) includes: Performing the following operations on each scattering point: (41), Assigning the shot - receiver pairs corresponding to the scattering point to different computing cores on the same node, and performing distributed parallel computing to obtain the local plane wave corresponding to each shot - receiver pair; (42), Superposing the local plane waves corresponding to all shot - receiver pairs to obtain the local point spread function corresponding to the scattering point.

6. The method for obtaining a point spread function with high cross - node parallelism according to claim 5, characterized in that: the operation of performing distributed parallel computing to obtain the local plane wave corresponding to each shot - receiver pair in step (41) includes: In each computing core, ray tracing is performed starting from the shot point to obtain the travel-time field T of the shot point s and the amplitude field A s , and ray tracing is performed starting from the geophone point to obtain the travel-time field T of the geophone point r and the amplitude field A r ; Calculating the angle using the following formula: where P S and P R are the incident directions of the shot point and the geophone respectively, θ is the angle, and x and z are the spatial X and Z axis directions; Calculating the local plane wave using the following formula: f′ = A s A r ∫ ||F(w)|| 2 exp(-iw(T s + T r ) + θ)dw where f′ is the local plane wave, w is the frequency, and F(w) is the result after Fourier transform of the input seismic wavelet.

7. The method for obtaining a point spread function with high cross - node parallelism according to claim 2, characterized in that: the operation of step (5) includes: Superposing the local point spread functions corresponding to all scattering points to obtain the global point spread function.

8. A system for obtaining a point spread function with high cross - node parallelism, characterized in that, The system comprises: An input unit, used to input the observation system, background velocity and seismic wavelet; A scattering point determination unit, connected to the input unit, for determining each scattering point in the point spread function model according to the reflection coefficient model; an allocation unit, connected to the scattering point determination unit, for allocating different scattering points in the point spread function model to different nodes; A local point spread function acquisition unit is connected to the allocation unit and is used to obtain a local point spread function corresponding to each scattering point; A global point spread function acquisition unit is connected to the local point spread function acquisition unit and is used to obtain a global point spread function; The output unit is connected to the global point spread function acquisition unit and is used to output the global point spread function.

9. A computer-readable storage medium, It is characterized in that The computer-readable storage medium stores at least one computer-executable program, and when the at least one program is executed by the computer, the computer executes the steps in the cross-node highly parallel point spread function acquisition method as described in any one of claims 1 to 7.

10. A computer device, It is characterized in that It comprises a memory and a processor, wherein the memory stores a computer program, and when the computer program is executed by the processor, the processor executes the steps of the cross-node highly parallel point spread function acquisition method as described in any one of claims 1 to 7.

Images