A method, system, device, and medium for multi-component seismic data migration imaging

Abstract

This invention discloses a method, system, device, and medium for multi-component seismic data migration imaging, relating to the fields of numerical simulation and imaging of seismic wave propagation. The method includes: discretizing a regular grid depth domain model of the seismic field; adaptively moving the initial discrete node position coordinates; constructing a finite difference numerical solution scheme for the elastic wave equation based on the final discrete node position coordinates; solving the finite difference numerical solution scheme for the elastic wave equation based on the final discrete node position coordinates; determining the pure P-wave field vector and pure S-wave field vector of the discrete nodes using a finite difference numerical calculation scheme for radial basis functions based on the source wavefield and detector wavefield; and performing imaging based on the pure P-wave field vector and pure S-wave field vector of the discrete nodes, according to the time consistency principle. This invention improves the accuracy and computational efficiency of migration imaging.

Description

Technical Field

[0001] This invention relates to the field of numerical simulation and imaging of seismic wave propagation, and in particular to a method, system, device and medium for multi-component seismic data migration imaging. Background Technology

[0002] Undulating terrain conditions are the first pressing challenge in seismic exploration. Dramatic topographic relief directly leads to a series of problems, including difficult data acquisition, low signal-to-noise ratio, severe scattering noise, inaccurate static correction, and low imaging accuracy. To eliminate the influence of undulating terrain, the most common practice in industry is to preprocess seismic data, such as performing static correction, before migration imaging. However, due to errors in these preprocessing methods, this process distorts the kinematic and dynamic characteristics of the seismic data, thus affecting the accuracy of subsequent migration imaging. Direct migration imaging based on undulating terrain is the fundamental method for solving seismic imaging problems in complex near-surface areas, and achieving efficient and high-precision numerical simulation of seismic waves under undulating terrain conditions is its core content. Although current research on numerical simulation of seismic waves under undulating terrain has yielded considerable results, a contradiction between accuracy and efficiency still exists in practical applications. Therefore, it is necessary to develop a set of efficient and high-precision theoretical methods for numerical simulation of seismic waves under undulating terrain to lay a solid wave theory foundation for direct migration imaging on undulating terrain.

[0003] Besides undulating surfaces, high-precision numerical simulation of seismic waves for complex structures is also crucial for seismic migration imaging. The classical finite difference method based on regular grids is currently the most widely used numerical simulation method for complex media in industry due to its simplicity and high efficiency. However, regular grids can cause the interface between the actual and numerical media parameters to not coincide, thus affecting the accuracy of the numerical simulation. Radial basis function (RBF) finite difference numerical simulation is a numerical solution algorithm for wave equations based on discrete, gridless nodes. This method is based on the radial basis function approximation theory, expressing the partial derivatives of the wave field as a linear combination of wave field values ​​from adjacent nodes. In contrast, the RBF method overcomes the grid dependency of the classical finite difference method, making it more suitable for solving seismic wave numerical simulation problems under complex structural conditions such as irregular geological bodies and undulating interfaces. However, influenced by factors such as node distribution, basis function type, deformation parameters, and boundary conditions, existing RBF methods face key challenges such as poor stability and low efficiency, making them difficult to apply to practical imaging processing. Therefore, researching a stable and efficient RBF method suitable for undulating surfaces and complex structures is of great significance for high-precision seismic imaging.

[0004] In conclusion, given the current complex geological conditions of both complex surface and complex structure, it is urgent to develop a migration imaging method that can accurately process seismic wave propagation under undulating surface conditions and achieve high efficiency and high precision for multi-component seismic data. Summary of the Invention

[0005] The purpose of this invention is to provide a method, system, device, and medium for multi-component seismic data migration imaging, which can improve the accuracy and computational efficiency of migration imaging.

[0006] To achieve the above objectives, the present invention provides the following solution:

[0007] A multi-component seismic data migration imaging method includes:

[0008] Obtain geological and geophysical data and surface elevation of the work area;

[0009] Based on the geological and geophysical data of the work area and the surface elevation of the work area, the regular grid depth domain model of the work area is discretized into nodes to obtain the initial discrete node position coordinates corresponding to the regular grid depth domain model.

[0010] The initial discrete node position coordinates are adaptively moved using the discrete node control equations to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model.

[0011] Based on the final discrete node position coordinates, a discrete node finite difference numerical solution scheme for the elastic wave equation is constructed using the radial basis function finite difference method.

[0012] The source wave field and detector wave field for each shot at each time step are obtained by using a filter to solve the discrete node finite difference numerical solution scheme of the elastic wave equation based on the final discrete node position coordinates.

[0013] Based on the source end wavefield and the detector wavefield, the pure P-wave field vector and pure S-wave field vector of the discrete node are determined using the discrete node elastic wavefield P-wave and S-wave separation radial basis function finite difference numerical calculation format.

[0014] Based on the principle of time consistency, the pure P-wave field vector and pure S-wave field vector of the discrete nodes are used to image the pure wave field vector at the source end and the pure wave field vector at the detector end using migration imaging conditions, so as to obtain the single-shot migration profile of each shot.

[0015] Optionally, based on the geological-geophysical data and surface elevation of the work area, the regular grid depth domain model of the work area is discretized into discrete nodes to obtain the initial discrete node position coordinates corresponding to the regular grid depth domain model, specifically including:

[0016] Using the geological and geophysical data of the work area and the physical boundary of the regular grid depth domain model of the work area based on the surface elevation of the work area, discrete nodes are subdivided layer by layer to obtain the initial spatial coordinates of the discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity and density at the discrete nodes.

[0017] Based on the initial spatial coordinates of the discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity, and density at the discrete nodes, the initial discrete node position coordinates corresponding to the regular mesh depth domain model are obtained using an adaptive meshless partitioning method.

[0018] Optionally, the initial discrete node position coordinates are adaptively moved using the discrete node control equations to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model, specifically including:

[0019] The initial discrete node position coordinates are solved using the discrete node control equations and the Gauss-Seidel iterative method to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model.

[0020] Optionally, the expression for the finite-difference numerical calculation scheme of the radial basis functions for separating the P-wave and S-waves in the discrete nodal elastic wave field is as follows:

[0021]

[0022]

[0023] in, v is the pure longitudinal wave field vector at discrete nodes. p The longitudinal wave velocity corresponding to the final discrete node position. The coefficients of the horizontal second derivative, u j For the horizontal component of the elastic wave field, For the finite difference coefficients of the radial basis functions with mixed second derivatives, w j This represents the vertical component of the elastic wave field. Let x be the unit vector in the x-direction of the Cartesian coordinate system. These are the coefficients of the vertical second derivative. Let be the unit vector in the z-direction of the Cartesian coordinate system. Let v be the pure transverse wave field vector at discrete nodes. s The transverse wave velocity corresponds to the final discrete node position, where i and j are both index numbers, and n is the total number of nodes participating in the difference calculation.

[0024] The present invention also provides a multi-component seismic data migration imaging system, comprising:

[0025] The acquisition module is used to acquire geological and geophysical data and surface elevation of the work area.

[0026] The discrete node subdivision module is used to perform discrete node subdivision on the regular grid depth domain model of the work area based on the geological-geophysical data of the work area and the surface elevation of the work area, so as to obtain the initial discrete node position coordinates corresponding to the regular grid depth domain model.

[0027] The node position adaptive movement module is used to adaptively move the node position based on the initial discrete node position coordinates using discrete node control equations, so as to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model.

[0028] The construction module is used to construct the discrete node finite difference numerical solution format of the elastic wave equation using the radial basis function finite difference method based on the final discrete node position coordinates;

[0029] The solution module is used to solve the discrete node finite difference numerical solution scheme of the elastic wave equation using a filter based on the final discrete node position coordinates, so as to obtain the source end wavefield and detector wavefield for each shot at each time.

[0030] The pure vector field determination module is used to determine the pure P-wave field vector and pure S-wave field vector of the discrete node based on the source end wave field and the detector wave field using the discrete node elastic wave field P-wave and S-wave separation radial basis function finite difference numerical calculation format.

[0031] The imaging module is used to image the pure wave field vector at the source end and the pure wave field vector at the detector end based on the principle of time consistency and the pure P-wave field vector and pure S-wave field vector of the discrete node using migration imaging conditions, so as to obtain the single-shot migration profile of each shot.

[0032] Optionally, the node position adaptive movement module specifically includes:

[0033] The layer-by-layer subdivision unit is used to perform discrete node subdivision layer by layer based on the physical boundary of the regular grid depth domain model of the work area, using the geological-geophysical data of the work area and the surface elevation of the work area based on the physical boundary of the work area, to obtain the initial spatial coordinates of the discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity and density at the discrete nodes.

[0034] An adaptive meshless partitioning element is used to obtain the initial discrete node position coordinates corresponding to the regular mesh depth domain model based on the initial spatial position coordinates of discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity and density at the discrete nodes, and the adaptive meshless partitioning method.

[0035] Optionally, the node position adaptive movement module specifically includes:

[0036] The node position adaptive moving unit is used to solve the initial discrete node position coordinates using the discrete node control equations and the Gauss-Seidel iterative method to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model.

[0037] Optionally, the expression for the finite-difference numerical calculation scheme of the radial basis functions for separating the P-wave and S-waves in the discrete nodal elastic wave field is as follows:

[0038]

[0039]

[0040] in, v is the pure longitudinal wave field vector at discrete nodes. p The longitudinal wave velocity corresponding to the final discrete node position. The coefficients of the horizontal second derivative, u j For the horizontal component of the elastic wave field, For the finite difference coefficients of the radial basis functions with mixed second derivatives, w j This represents the vertical component of the elastic wave field. Let x be the unit vector in the x-direction of the Cartesian coordinate system. These are the coefficients of the vertical second derivative. Let be the unit vector in the z-direction of the Cartesian coordinate system. Let v be the pure transverse wave field vector at discrete nodes. s The transverse wave velocity corresponds to the final discrete node position, where i and j are both index numbers, and n is the total number of nodes participating in the difference calculation.

[0041] The present invention also provides an electronic device, comprising:

[0042] One or more processors;

[0043] A storage device on which one or more programs are stored;

[0044] When the one or more programs are executed by the one or more processors, the one or more processors cause the one or more processors to perform the method as described.

[0045] The present invention also provides a computer storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the method as described.

[0046] According to specific embodiments provided by the present invention, the present invention discloses the following technical effects:

[0047] This invention acquires geological and geophysical data and surface elevation of the work area; based on the geological and geophysical data and surface elevation, it performs discrete node subdivision on a regular grid depth domain model of the work area to obtain the initial discrete node position coordinates corresponding to the regular grid depth domain model; it then uses discrete node control equations to adaptively move the node positions based on the initial discrete node position coordinates to obtain the final discrete node position coordinates corresponding to the regular grid depth domain model; based on the final discrete node position coordinates, it constructs a discrete node finite difference numerical solution format for the elastic wave equation using the radial basis function finite difference method; and finally, it uses a filter to solve the problem based on the final discrete node position coordinates. The discrete node position coordinates are used to solve the elastic wave equation using the discrete node finite difference numerical solution scheme to obtain the source end wavefield and detector wavefield for each shot at each time step. Based on the source end wavefield and detector wavefield, the pure P-wave field vector and pure S-wave field vector of the discrete node are determined using the discrete node elastic wave field P-wave and S-wave separation radial basis function finite difference numerical calculation scheme. Based on the time consistency principle, the pure P-wave field vector at the source end and the pure S-wave field vector at the detector end are imaged using migration imaging conditions based on the pure P-wave field vector and pure S-wave field vector of the discrete node to obtain the single-shot migration profile for each shot, thereby improving the accuracy and computational efficiency of migration imaging. Attached Figure Description

[0048] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0049] Figure 1 Flowchart of the multi-component seismic data migration imaging method provided by the present invention;

[0050] Figure 2 A model diagram of a uniform medium with undulating terrain;

[0051] Figure 3 for Figure 2 The diagram shows the adaptive undulating surface discrete node partitioning result corresponding to the uniform medium undulating surface model.

[0052] Figure 4 A 0.4s seismic wave field snapshot of a convex Gaussian undulating surface model;

[0053] Figure 5 A 0.4s seismic wave field snapshot of a concave Gaussian undulating surface model;

[0054] Figure 6 To utilize Figure 4The image shows a seismic record simulated by a convex Gaussian undulation model.

[0055] Figure 7 To utilize Figure 5 The image shows a seismic record simulated by a concave Gaussian undulation model.

[0056] Figure 8 This is a model diagram of an undulating surface depression.

[0057] Figure 9 for Figure 8 The discrete node partitioning result of the undulating surface depression model is shown in the figure.

[0058] Figure 10 In order to be in Figure 8 A snapshot of the seismic wavefield obtained using the method of this invention in the undulating surface depression model shown;

[0059] Figure 11 In order to be in Figure 8 The seismic record map obtained using the method of this invention is shown in the undulating surface depression model.

[0060] Figure 12 for Figure 8 The image shows a multi-shot migration imaging profile of the undulating surface depression model.

[0061] Figure 13 A model diagram of the undulating Marmousi-2 terrain;

[0062] Figure 14 In order to be in Figure 13 The image shown is a snapshot of the multi-component seismic wavefield obtained using the method of this invention in the Marmousi-2 model of the undulating terrain.

[0063] Figure 15 In order to be in Figure 13 A snapshot of the pure P- and S-wave fields obtained using the method of this invention in the Marmousi-2 model of the undulating terrain shown.

[0064] Figure 16 In order to be in Figure 13 The image shown is a multi-shot migration imaging profile obtained using the method of this invention in the Marmousi-2 model of undulating terrain. Detailed Implementation

[0065] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0066] The purpose of this invention is to provide a method, system, device, and medium for multi-component seismic data migration imaging, which can improve the accuracy and computational efficiency of migration imaging.

[0067] To make the above-mentioned objects, features and advantages of the present invention more apparent and understandable, the present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0068] This addresses the key limitation of current finite difference methods, which are ill-suited for high-precision numerical simulations and migration imaging in complex dual-medium media. For example... Figure 1 As shown, the present invention provides a multi-component seismic data migration imaging method, comprising:

[0069] Step 101: Obtain geological and geophysical data and surface elevation of the work area.

[0070] Step 102: Based on the geological and geophysical data of the work area and the surface elevation of the work area, the regular grid depth domain model of the work area is discretized into nodes to obtain the initial discrete node position coordinates corresponding to the regular grid depth domain model.

[0071] Step 102 specifically includes: using the geological-geophysical data of the work area and the surface elevation of the work area based on the physical boundary of the regular grid depth domain model of the work area, performing discrete node subdivision layer by layer to obtain the initial spatial coordinates of the discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity and density at the discrete nodes; based on the initial spatial coordinates of the discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity and density at the discrete nodes, and the adaptive meshless subdivision method, obtaining the initial discrete node position coordinates corresponding to the regular grid depth domain model.

[0072] Based on the geological and geophysical characteristics of the work area, a regular grid depth domain model of the work area is established. This model specifically includes a velocity model and a density model ρ(x). At the same time, the observation system parameters, source parameters, and wavefield extension parameters for migration imaging are determined. Among them, the velocity model and density model are existing models. The only task is to extract the existing models and obtain the depth domain model based on the specific geological and geophysical characteristics of the work area.

[0073] The velocity model includes the longitudinal wave velocity v. p (x) and transverse wave velocity v s (x), where x = (x, z) represents the spatial position vector in the Cartesian coordinate system, and x and z represent its horizontal and vertical components.

[0074] The specific parameters of the observation system include: the number of seismic sources, the location of each seismic source and the corresponding spatial distribution characteristics of the geophones, the spatial distribution characteristics of the seismic sources, the geophone spacing, the maximum offset distance, and the minimum offset distance.

[0075] The source parameters specifically include: source wavelet type, source wavelet dominant frequency, and source spatial characteristics.

[0076] The wavefield extension parameters specifically include: time step, time sampling length, absorbing boundary layer thickness, and wavefield extension accuracy.

[0077] The surface elevation of the work area is obtained. Using the surface elevation information and the physical space range of the model, the velocity model and density model ρ(x) of the regular grid depth domain model are discretized into discrete nodes to obtain the initial discrete node position coordinates corresponding to the depth domain model.

[0078] The regular mesh deep domain model implements discrete node partitioning, specifically including:

[0079] (a) Based on the physical boundary of the depth domain model, starting from the lower left boundary position of the model, discrete nodes are decomposed layer by layer in a clockwise direction. The node radius r of adjacent nodes is determined by the transverse wave velocity and longitudinal wave velocity at the corresponding physical position. Finally, the initial spatial position coordinates of the discrete nodes on the physical boundary, as well as the longitudinal wave velocity, transverse wave velocity and density at the nodes are obtained.

[0080] (b) Based on the adaptive meshless partitioning method, starting from the initial discrete nodes on the physical boundary described in step (a), the node radius d of the adjacent discrete nodes in the next layer is determined based on the transverse wave velocity and longitudinal wave velocity at the corresponding physical location. This process is repeated layer by layer until the model space is fully covered. Combining the discrete nodes on the physical boundary described in step (a), the initial spatial coordinates x0 = (η, ξ) of the discrete nodes corresponding to the depth domain model and the longitudinal wave velocity v at the nodes are finally obtained. p (x0), transverse wave velocity v s (x0) and density ρ(x0), where η and ξ represent the horizontal and vertical coordinate components, respectively, in the Cartesian coordinate system, which are different from x and z.

[0081] (c) The result of step (b) is the node radius of each discrete node relative to its neighboring nodes. Nodes with radii smaller than d are eliminated. min / 2 nodes, to obtain the initial spatial position coordinates x0 of discrete nodes with a more uniform distribution;

[0082] d min This represents the minimum value of the discrete node radius corresponding to the depth domain model, which is determined by the following formula:

[0083]

[0084] In formula (1), f represents the minimum shear wave velocity value corresponding to the depth domain model. max This indicates the maximum frequency of the wavelet.

[0085] Step 103: Adaptively move the node positions using the discrete node control equations on the initial discrete node position coordinates to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model.

[0086] Step 103 specifically includes: solving the initial discrete node position coordinates using the discrete node control equations and the Gauss-Seidel iterative method to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model.

[0087] Based on the initial discrete node position coordinates, the node positions are adaptively moved using the discrete node control equations to obtain the final discrete node position coordinates X = (X, Z) corresponding to the depth domain model. Then, inverse distance interpolation is performed based on the velocity and density of the original regular mesh depth domain model.

[0088] Simultaneously, the longitudinal wave velocity v corresponding to the final discrete node position is obtained. p (X), transverse wave velocity v s (X) and density ρ(X).

[0089] The specific governing equations for discrete nodes are as follows:

[0090]

[0091] In equation (2), X and Z represent the horizontal and vertical components of the position coordinate X, and the weighting function ω is as follows:

[0092]

[0093] Where ε is an empirical parameter of the weighting function.

[0094] In equation (2), the intermediate parameters a1, a2, a3, b1, b2, b3, c1, c2, c3 are as follows:

[0095]

[0096] ξ and η are the components of the deformation of X and Z in the x and z directions of the Cartesian coordinate system, respectively. In equation (2), J represents the Jacobi polynomial, specifically:

[0097]

[0098] Among them, the adaptive movement of node positions is implemented using discrete node control equations, with the following specific characteristics:

[0099] Based on the aforementioned discrete node control equation, the Gauss-Seidel iterative method is used to solve the equation, thereby obtaining the optimal solution of the discrete node control equation, which is the final position coordinate X of the discrete node corresponding to the depth domain model.

[0100] Step 104: Based on the final discrete node position coordinates, construct the discrete node finite difference numerical solution format for the elastic wave equation using the radial basis function finite difference method.

[0101] Based on the final discrete node position coordinates, a discrete node finite difference numerical solution scheme for the elastic wave equation is constructed using the radial basis function finite difference method, as follows:

[0102]

[0103] In equation (6), u and w represent the horizontal and vertical components of the elastic wave field, respectively; subscripts i and j represent the index numbers of the current central node and the difference node used for numerical solution; superscript l represents time; τ represents the time sampling step; n represents the total number of nodes participating in the difference calculation; λ(X) and μ(X) represent the Lamé constants corresponding to the discrete node model; and subscript i also represents the retrieval number of the current central node. and denoted by finite difference coefficients of radial basis functions for the horizontal second derivative, vertical second derivative, and mixed second derivative, respectively.

[0104] in It is obtained by solving the following system of equations, specifically:

[0105]

[0106] It is obtained by solving the following system of equations, specifically:

[0107]

[0108] It is obtained by solving the following system of equations, specifically:

[0109]

[0110] In equations (7), (8) and (9), X i This represents the vector representing the current center node position. The radial basis functions are:

[0111]

[0112] In equation (10), δ represents the radial basis function. The deformation parameter is usually taken as 0.02; the radial basis function independent variable r, which is the discrete node radius, is specifically r = ||X i -X j || represents the distance between the central node with index i and its neighboring differential node with index j.

[0113] Step 105: Using a filter, solve the discrete node finite difference numerical solution scheme of the elastic wave equation according to the final discrete node position coordinates to obtain the source end wavefield and detector wavefield for each shot at each time.

[0114] For each shot and each time step, perform the following operations: For each shot, at each time step of that shot, obtain the source location, source wavelet, and source spatial characteristic function of that shot based on the observation system parameters, obtain the corresponding source function by multiplying the two, and set the source wavelet at the corresponding source location.

[0115] Based on the P-wave velocity v corresponding to the final discrete node position p (X), transverse wave velocity v s (X) and density ρ(X), the above elastic wave equation is solved by the discrete node finite difference numerical formula (6), the elastic wave equation is solved, the positive extension of the source is realized, and the source end wave field of the source at each moment is obtained.

[0116] Similarly, for each shot, at each time step of that shot, the spatial distribution of the detector array corresponding to that shot is obtained based on the observation system parameters. Along the counter-clockwise direction, the multi-component seismic data D(x,t) of each time step of that shot is used as boundary conditions, based on the P-wave velocity v corresponding to the final discrete node position. p (X), transverse wave velocity v s (X) and density ρ(X), the elastic wave equation is solved by the discrete node finite difference numerical formula (6) of the above elastic wave equation, and the corresponding detector wave field is reversed to obtain the detector end wave field at each moment of the source, and finally the source end wave field and detector end wave field at each moment of each shot are obtained.

[0117] In each time step of the forward extension of the wave field at the source end and the reverse extension of the wave field at the detector end, the instability of the wave field is suppressed by applying a stability filter to the extended wave field. The application of the stability filter is to first apply the filtering operation shown in formula (11) to the wave field at the current time l in each time step, and then obtain the new filtered wave field at the current time. Then, this new filtered wave field replaces the wave field at the current time in the original formula (6).

[0118] Specifically, the stability filter is characterized by performing the following operation on the extended wavefield at each time step:

[0119]

[0120] In equation (11), β x and β z This represents the filter parameter, typically set to 0.02. and Let G and H represent the horizontal and vertical components of the elastic wave field after processing by the stability filter, respectively, with intermediate variables G and H being:

[0121]

[0122] In each time step of the forward extension of the wavefield at the source end and the reverse extension of the wavefield at the receiver end, a stability filter is applied to the extended wavefield to suppress the instability artifacts of the wavefield. Specifically, the horizontal and vertical components of the elastic wavefield processed by the stability filter are used to replace the original horizontal and vertical components of the elastic wavefield for wavefield extension. That is, by substituting formula (12) into formula (11), the following can be obtained. According to Substituting into formula (6) yields formula (13), which can be specifically expressed as:

[0123]

[0124] Step 106: Based on the source end wavefield and the detector wavefield, determine the pure P-wave field vector and pure S-wave field vector of the discrete node using the discrete node elastic wavefield P-wave separation radial basis function finite difference numerical calculation format.

[0125] For each shot at each time step, perform the following operations: At each time step of each shot, based on the accurate P-wave and S-wave separation equations, construct a finite-difference numerical calculation scheme for the radial basis functions of the discrete node elastic wavefield P-wave and S-wave separation. Apply this numerical calculation scheme to the seismic wavefields at the source and detector ends obtained in step 105 to obtain the corresponding pure P-wave field vectors of the discrete nodes. and pure transverse wave field vector

[0126] The precise equations for separating the P-waves and S-waves are as follows:

[0127]

[0128] Among them, U P =(u P ,w P ) TU represents the precise pure longitudinal wave field vector. S =(u S ,w S ) T The precise pure transverse wave field vector is represented by the scalar longitudinal wave field P and the vector transverse wave field S in formula (14), which are specifically:

[0129]

[0130] Among them, U=(u,w) T This represents the seismic wave field vector.

[0131] Using the radial basis function finite difference method, equation (14) is numerically discretized to construct a numerical calculation scheme for the radial basis function finite difference method for separating the P-wave and S-wave of the discrete nodal elastic wave field, specifically as follows:

[0132]

[0133] in, and These represent the unit direction vectors in the x and z directions, respectively, in the Cartesian coordinate system. Let be the pure longitudinal wave field vector of the i-th discrete node at time l. Let be the pure transverse wave field vector of the i-th discrete node at time l. The finite difference coefficients of the radial basis functions with horizontal second derivatives are... The finite difference coefficients of the radial basis functions with the vertical second derivative are... The finite difference coefficients are the radial basis functions of the mixed second derivative.

[0134] Operators in formulas (14) and (15) and These represent gradient, divergence, and curl operations, respectively, and contain different spatial partial derivatives. The spatial partial derivatives are numerically represented by the radial basis function finite difference, thus obtaining formula (16). In formula (16), u and w represent the horizontal and vertical components of the elastic wave field, respectively. Both the source end wave field and the wave generator end wave field contain these two components. By applying formula (16) respectively, the aforementioned purpose can be achieved.

[0135] Step 107: Based on the principle of time consistency, the pure P-wave field vector and pure S-wave field vector of the discrete node are used to image the pure wave field vector at the source end and the pure wave field vector at the detector end using the migration imaging conditions to obtain the single-shot migration profile of each shot.

[0136] At the same time, based on the principle of time consistency, the migration imaging condition is applied to the pure wavefield vector obtained in step 106, wherein the pure wavefield vector includes the pure P-wave field vector and the pure S-wave field vector; the pure wavefield vector s at the source end is obtained for each time step of each shot. m and the pure wave field vector r at the detector end n Imaging is performed to obtain the single-shot offset profile of each shot; then the single-shot offset profiles of all shots are obtained. Based on the coverage relationship contained in the observation system, the single-shot offset profiles of all shots are superimposed to obtain the final multi-component offset profile.

[0137] Specifically, the offset imaging conditions are as follows:

[0138] I mn (x)=∫s m (x,t)·r n (x,t)dt,m,n=p,s (17)

[0139] Among them, I mn Indicates the offset profile, s m and r n Let represent the pure wave field vector m-wave mode at the source end and the pure wave field vector n-wave mode at the detector end, respectively, where t is time.

[0140] The pure wave field vector s at the source end m Specifically including s P and s S The pure wave field vector r at the detector end n Specifically, r P and r S The terms s and r represent the seismic wave field vector at the source end and the seismic wave field vector at the detector end, respectively.

[0141] The present invention, by adopting the above technical solutions, has the following advantages: 1) The method of the present invention is a migration imaging method for multi-component seismic data in dual complex media. Compared with the traditional migration imaging method based on the conventional finite difference method, the method of the present invention can accurately handle undulating surface conditions, and the spacing between discrete nodes automatically changes with the parameter model, resulting in higher computational efficiency and accuracy; 2) The method of the present invention provides a numerical calculation method for wave equations based on discrete nodes. Compared with existing methods, the method of the present invention does not require manual selection of deformation parameters, and effectively suppresses the instability phenomenon in the wave field simulation process based on discrete nodes through a stability filter; 3) Compared with conventional curved coordinate system algorithms, the difference form of the present invention is simpler, programming is easier, parallel algorithm design is simpler, and it is easier to promote and apply in production; 4) The present invention provides a good technical development direction for migration imaging and waveform inversion algorithms for complex structures of undulating surfaces in the future, and has great application potential.

[0142] Figure 2 It is a model of a undulating surface in a uniform medium: where, Figure 2 (a) in the model is a convex Gaussian undulation model. Figure 2 (b) in the model is a concave Gaussian undulation model. Figure 3 yes Figure 2 The adaptive undulating surface discrete node partitioning result corresponding to the uniform medium undulating surface model shown is as follows: Figure 3 In the diagram, (a) shows the node partitioning results of the convex Gaussian undulation model. Figure 3 (b) in the figure is the node subdivision result of the concave Gaussian undulation model. It can be seen that both exhibit the characteristic of the same node spacing. Figure 4 This is a 0.4s seismic wave field snapshot of a convex Gaussian undulating surface model: where, Figure 4 In the figure, (a) represents the horizontal component obtained by the method of the present invention. Figure 4 (b) in the figure represents the vertical component obtained by the method of the present invention. Figure 4 In the diagram, (c) represents the level component obtained using the conventional finite difference method. Figure 4 In this context, (d) represents the vertical component obtained using the conventional finite difference method. Figure 4 As can be seen, compared with conventional finite difference methods, the method of the present invention does not produce step scattering in the undulating structure, the simulation results are more accurate, and the phase axis of the reflected wave is clearer. Figure 5 This is a 0.4s seismic wave field snapshot of a concave Gaussian undulating surface model: where, Figure 5 In the figure, (a) represents the horizontal component obtained by the method of the present invention. Figure 5 (b) in the figure represents the vertical component obtained by the method of the present invention. Figure 5 In the diagram, (c) represents the level component obtained using the conventional finite difference method. Figure 5 In the figure, (d) represents the vertical component obtained by the conventional finite difference method, derived from... Figure 5 Similar conclusions can be drawn from the convex Gaussian undulation surface model. Figure 6 It is to utilize Figure 4 The earthquake record simulated by the convex Gaussian undulation model shown is as follows: Figure 6 (a) in the figure represents the horizontal component obtained using the method of this invention. Figure 6 (b) The vertical component obtained by the method of the present invention. Figure 6 (c) in the figure represents the level component obtained using the conventional finite difference method. Figure 6 (d) is the level component obtained using the traditional finite difference method. Figure 7 It is to utilize Figure 5 The earthquake record simulated by the concave Gaussian undulation model shown is as follows: Figure 7 (a) in the figure represents the horizontal component obtained using the method of this invention. Figure 7 (b) The vertical component obtained by the method of the present invention. Figure 7 (c) in the figure represents the level component obtained using the conventional finite difference method. Figure 7 (d) shows the horizontal component obtained using the traditional finite difference method. It is clear from this that the method of the present invention has advantages over the conventional finite difference method in eliminating staircase scattering, and the surface wave phase axis is clear. The obtained wave field is very stable at the undulating interface and does not produce a large amount of secondary scattering. It can accurately perform numerical simulation of undulating surface models under free boundary conditions.

[0143] Figure 8 It is an undulating surface depression model: where, Figure 8 In the figure, (a) represents the longitudinal wave velocity. Figure 8 In the equation (b), the transverse wave velocity is represented. Figure 9 yes Figure 8 The discrete node partitioning results of the undulating surface depression model shown are as follows: Figure 9 (a) in the figure is the node subdivision result of the method of the present invention. Figure 9 (b) in the figure is the result of conventional node partitioning. Compared with conventional node partitioning methods, the node spacing in the node partitioning method of this invention increases with the increase of speed. Figure 10 Is Figure 8 A snapshot of the seismic wavefield obtained using the method of this invention in the undulating surface depression model shown: where, Figure 10 (a) in the figure represents the 0.5s horizontal component. Figure 10 (b) in the figure represents the 0.5s vertical component. Figure 10 (c) in the figure represents the horizontal component at 1.0 s. Figure 10 In this context, (d) represents the vertical component at 1.0s. Figure 10 (e) in the equation represents the 1.5s horizontal component. Figure 10The (f) 1.5s vertical component. Comparison of seismic wave field snapshots at different times shows that the method of this invention does not exhibit step scattering in the near-surface region, and the wave remains stable and free of dispersion when propagating to deeper regions. Figure 11 Is Figure 8 The seismic record obtained using the method of this invention in the undulating surface depression model shown is as follows: Figure 11 (a) in the text represents the horizontal component. Figure 11 (b) in the diagram represents the vertical component. (From...) Figure 11 As can be seen, the method of the present invention can handle undulating surfaces well, and the direct wave, reflected wave, and converted wave can be clearly seen, which verifies the effectiveness of the method of the present invention. Figure 12 yes Figure 8 The multi-shot migration imaging profile of the undulating surface depression model shown is as follows: Figure 12 (a) in the image is the PP imaging profile. Figure 12 (b) in the figure is a PS imaging profile, which shows that the method of the present invention can flexibly handle undulating surfaces, and the location of the subsurface interface is accurate, with clear and continuous phase axes. The numerical results verify the correctness and effectiveness of the method of the present invention.

[0144] Figure 13 It is the Marmousi-2 model of undulating terrain: where, Figure 13 In the figure, (a) represents the longitudinal wave velocity. Figure 13 (b) in the figure represents the transverse wave velocity. Numerical simulations were performed based on this model with a time step of 0.5 ms and a recording time of 4.5 s. The Ricker wavelet with a dominant frequency of 20 Hz was selected as the source time function. Detector points were uniformly distributed on the ground surface with a receiving range of 10 km and a spacing of 10 m between adjacent detector points. Figure 14 Is Figure 13 The snapshot of the multi-component seismic wavefield obtained using the method of this invention in the Marmousi-2 model showing an undulating surface is shown below: Figure 14 (a) in the figure represents the horizontal component at 1.0 s. Figure 14 (b) in the figure represents the vertical component at 1.0s. Figure 14 (c) in the figure represents the 2.0s horizontal component. Figure 14 (d) in the figure represents the 2.0s vertical component. Figure 14 The results show that the method of the present invention can achieve high-precision numerical simulation for complex surface structures with drastic changes in surface elevation. Figure 15 Is Figure 13 The snapshot of the pure P- and S-wave fields obtained using the method of this invention in the Marmousi-2 model of the undulating surface shown is as follows: Figure 15 (a) in the image is a 1.0s pure P-wave snapshot. Figure 15 (b) in the image is a 1.0s pure shear wave snapshot. Figure 15(c) in the image is a 2.0s pure P-wave snapshot. Figure 15 (d) in the figure is a 2.0s pure transverse wave snapshot. As can be seen from the figure, the method of the present invention can effectively separate pure longitudinal and transverse wave fields. Figure 16 Is Figure 13 The multi-shot migration imaging profile obtained using the method of this invention in the Marmousi-2 model showing undulating terrain is shown below: where, Figure 16 (a) in the image is the PP imaging profile. Figure 16 (b) in the figure is the PS imaging profile. It can be seen that the imaging results are consistent with the construction interface information of the migration velocity model, and the location of the subsurface interface is accurate. Whether it is PP imaging or PS imaging, the phase axis is clear and continuous, and the signal-to-noise ratio is high, indicating that the method of the present invention has good adaptability to dual complex medium models.

[0145] The present invention also provides a multi-component seismic data migration imaging system, comprising:

[0146] The acquisition module is used to acquire geological and geophysical data and surface elevation of the work area.

[0147] The discrete node subdivision module is used to perform discrete node subdivision on the regular grid depth domain model of the work area based on the geological-geophysical data and the surface elevation of the work area, so as to obtain the initial discrete node position coordinates corresponding to the regular grid depth domain model.

[0148] The node position adaptive movement module is used to adaptively move the node position based on the initial discrete node position coordinates using the discrete node control equations, so as to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model.

[0149] The module is used to construct a discrete node finite difference numerical solution format for the elastic wave equation using the radial basis function finite difference method based on the final discrete node position coordinates.

[0150] The solution module is used to solve the discrete node finite difference numerical solution scheme of the elastic wave equation using a filter based on the final discrete node position coordinates, so as to obtain the source end wavefield and detector wavefield for each shot at each time.

[0151] The pure vector field determination module is used to determine the pure P-wave field vector and pure S-wave field vector of the discrete node based on the source end wave field and the detector wave field using the discrete node elastic wave field P-wave and S-wave separation radial basis function finite difference numerical calculation format.

[0152] The imaging module is used to image the pure wave field vector at the source end and the pure wave field vector at the detector end based on the principle of time consistency and the pure P-wave field vector and pure S-wave field vector of the discrete node using migration imaging conditions, so as to obtain the single-shot migration profile of each shot.

[0153] As an optional implementation, the node position adaptive movement module specifically includes:

[0154] The layer-by-layer subdivision unit is used to perform discrete node subdivision layer by layer based on the physical boundary of the regular grid depth domain model of the work area, using the geological-geophysical data of the work area and the surface elevation of the work area. This process yields the initial spatial coordinates of the discrete nodes on the physical boundary, as well as the P-wave velocity, S-wave velocity, and density at the discrete nodes.

[0155] An adaptive meshless partitioning element is used to obtain the initial discrete node position coordinates corresponding to the regular mesh depth domain model based on the initial spatial position coordinates of discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity and density at the discrete nodes, and the adaptive meshless partitioning method.

[0156] As an optional implementation, the node position adaptive movement module specifically includes:

[0157] The node position adaptive moving unit is used to solve the initial discrete node position coordinates using the discrete node control equations and the Gauss-Seidel iterative method to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model.

[0158] As an optional implementation, the expression for the finite-difference numerical calculation scheme of the radial basis functions for separating the P-wave and S-waves in the discrete nodal elastic wave field is as follows:

[0159]

[0160]

[0161] in, v is the pure longitudinal wave field vector at discrete nodes. p The longitudinal wave velocity corresponding to the final discrete node position. The coefficients of the horizontal second derivative, u j For the horizontal component of the elastic wave field, For the finite difference coefficients of the radial basis functions with mixed second derivatives, w j This represents the vertical component of the elastic wave field. Let x be the unit vector in the x-direction of the Cartesian coordinate system. These are the coefficients of the vertical second derivative. Let be the unit vector in the z-direction of the Cartesian coordinate system. Let v be the pure transverse wave field vector at discrete nodes. s The transverse wave velocity corresponds to the final discrete node position, where i and j are both index numbers, and n is the total number of nodes participating in the difference calculation.

[0162] The present invention provides an electronic device comprising: one or more processors; a storage device having one or more programs stored thereon; wherein when the one or more programs are executed by the one or more processors, the one or more processors perform the method as described above.

[0163] The present invention also provides a computer storage medium having a computer program stored thereon, wherein the computer program, when executed by a processor, implements the method as described.

[0164] This invention relates to a migration imaging method for multi-component seismic data in dual complex media. The method includes: establishing a depth-domain velocity and density model based on the geological and geophysical characteristics of the target area; using actual surface elevation information to perform discrete node subdivision on the obtained model to obtain initial discrete nodes; using the discrete node control equations to perform adaptive migration and obtain the final discrete node coordinates; based on the final discrete nodes, constructing a discrete node finite difference numerical solution scheme for the elastic wave equation; based on a stability filter and given the source type, performing a stable forward extension of the elastic wave field at the source end of the discrete nodes, and simultaneously performing a stable reverse time extension of the elastic wave field at the detector end of the multi-component seismic data; based on an accurate P-wave and S-wave separation equation, constructing a P-wave and S-wave separation difference scheme for the discrete node elastic wave field, and applying it to the elastic wave fields at the source end and detector end to obtain the pure P-wave field and pure S-wave field of the discrete nodes; and using the obtained pure wave field to perform imaging to obtain a migration profile. This invention can effectively simulate the propagation laws and response characteristics of seismic waves in dual complex media and realize migration imaging of multi-component seismic data. This invention provides effective theoretical and technical support for revealing the propagation mechanism of deep / ultra-deep complex seismic waves under undulating surface conditions and for imaging multi-component seismic data.

[0165] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on its differences from other embodiments. Similar or identical parts between embodiments can be referred to interchangeably. For the systems disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the descriptions are relatively simple; relevant parts can be referred to the method section.

[0166] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A method for multi-component seismic data migration imaging, characterized in that, include: Obtain geological and geophysical data and surface elevation of the work area; Based on the geological and geophysical data of the work area and the surface elevation of the work area, the regular grid depth domain model of the work area is discretized into nodes to obtain the initial discrete node position coordinates corresponding to the regular grid depth domain model. The initial discrete node position coordinates are adaptively moved using the discrete node control equations to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model. Based on the final discrete node position coordinates, a discrete node finite difference numerical solution scheme for the elastic wave equation is constructed using the radial basis function finite difference method. The discrete-node finite-difference numerical solution scheme of the elastic wave equation is solved using a filter based on the final discrete node position coordinates to obtain the source wavefield and detector wavefield for each shot at each time step. The filter is a stability filter, specifically characterized by performing the following operations on the extended wavefield at each time step: ； and Indicates filter parameters, and These represent the horizontal and vertical components of the elastic wave field after processing by the stability filter, respectively. l For a specific moment; Indicates the time sampling step size. Indicates the index number of the current central node; Based on the source end wavefield and the detector wavefield, the pure P-wave field vector and pure S-wave field vector of the discrete node are determined using the discrete node elastic wavefield P-wave and S-wave separation radial basis function finite difference numerical calculation format. Based on the principle of time consistency, the pure P-wave field vector and pure S-wave field vector of the discrete nodes are used to image the pure wave field vector at the source end and the pure wave field vector at the detector end using migration imaging conditions, so as to obtain the single-shot migration profile of each shot.

2. The multi-component seismic data migration imaging method according to claim 1, characterized in that, Based on the geological and geophysical data and surface elevation of the work area, the regular grid depth domain model of the work area is discretized into discrete nodes to obtain the initial discrete node position coordinates corresponding to the regular grid depth domain model, specifically including: Using the geological and geophysical data of the work area and the physical boundary of the regular grid depth domain model of the work area based on the surface elevation of the work area, discrete nodes are subdivided layer by layer to obtain the initial spatial coordinates of the discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity and density at the discrete nodes. Based on the initial spatial coordinates of the discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity, and density at the discrete nodes, the initial discrete node position coordinates corresponding to the regular mesh depth domain model are obtained using an adaptive meshless partitioning method.

3. The multi-component seismic data migration imaging method according to claim 1, characterized in that, The initial discrete node position coordinates are adaptively moved using the discrete node control equations to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model. Specifically, this includes: The initial discrete node position coordinates are solved using the discrete node control equations and the Gauss-Seidel iterative method to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model.

4. The multi-component seismic data migration imaging method according to claim 1, characterized in that, The expression for the finite difference numerical calculation scheme of the radial basis function for separating the longitudinal and transverse waves in the discrete node elastic wave field is as follows: in, Let be the pure longitudinal wave field vector at the discrete nodes. The longitudinal wave velocity corresponding to the final discrete node position. These are the coefficients of the horizontal second derivative. For the horizontal component of the elastic wave field, The finite difference coefficients of the radial basis functions with mixed second derivatives. This represents the vertical component of the elastic wave field. In Cartesian coordinate system x The unit vector of direction, These are the coefficients of the vertical second derivative. In Cartesian coordinate system z The unit vector of direction, For the pure transverse wave field vector of the discrete nodes, The transverse wave velocity corresponding to the final discrete node position. i and j All are index numbers. n This represents the total number of nodes participating in the differential calculation.

5. A multi-component seismic data migration imaging system, characterized in that, include: The acquisition module is used to acquire geological and geophysical data and surface elevation of the work area. The discrete node subdivision module is used to perform discrete node subdivision on the regular grid depth domain model of the work area based on the geological-geophysical data of the work area and the surface elevation of the work area, so as to obtain the initial discrete node position coordinates corresponding to the regular grid depth domain model. The node position adaptive movement module is used to adaptively move the node position based on the initial discrete node position coordinates using discrete node control equations, so as to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model. The construction module is used to construct the discrete node finite difference numerical solution format of the elastic wave equation using the radial basis function finite difference method based on the final discrete node position coordinates; The solution module is used to solve the discrete-node finite-difference numerical solution scheme of the elastic wave equation using a filter based on the final discrete node position coordinates, to obtain the source wavefield and detector wavefield for each shot at each time step; the filter is a stability filter, specifically characterized by performing the following operations on the extended wavefield at each time step: ； and Indicates filter parameters, and These represent the horizontal and vertical components of the elastic wave field after processing by the stability filter, respectively. l For a specific moment; Indicates the time sampling step size. Indicates the index number of the current central node; The pure vector field determination module is used to determine the pure P-wave field vector and pure S-wave field vector of the discrete node based on the source end wave field and the detector wave field using the discrete node elastic wave field P-wave and S-wave separation radial basis function finite difference numerical calculation format. The imaging module is used to image the pure wave field vector at the source end and the pure wave field vector at the detector end based on the principle of time consistency and the pure P-wave field vector and pure S-wave field vector of the discrete node using migration imaging conditions, so as to obtain the single-shot migration profile of each shot.

6. The multi-component seismic data migration imaging system according to claim 5, characterized in that, The node position adaptive movement module specifically includes: The layer-by-layer subdivision unit is used to perform discrete node subdivision layer by layer based on the physical boundary of the regular grid depth domain model of the work area, using the geological-geophysical data of the work area and the surface elevation of the work area based on the physical boundary of the work area, to obtain the initial spatial coordinates of the discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity and density at the discrete nodes. An adaptive meshless partitioning element is used to obtain the initial discrete node position coordinates corresponding to the regular mesh depth domain model based on the initial spatial position coordinates of discrete nodes on the physical boundary, the P-wave velocity, S-wave velocity and density at the discrete nodes, and the adaptive meshless partitioning method.

7. The multi-component seismic data migration imaging system according to claim 5, characterized in that, The node position adaptive movement module specifically includes: The node position adaptive moving unit is used to solve the initial discrete node position coordinates using the discrete node control equations and the Gauss-Seidel iterative method to obtain the final discrete node position coordinates corresponding to the regular mesh depth domain model.

8. The multi-component seismic data migration imaging system according to claim 5, characterized in that, The expression for the finite difference numerical calculation scheme of the radial basis function for separating the longitudinal and transverse waves in the discrete node elastic wave field is as follows: in, Let be the pure longitudinal wave field vector at the discrete nodes. The longitudinal wave velocity corresponding to the final discrete node position. These are the coefficients of the horizontal second derivative. For the horizontal component of the elastic wave field, The finite difference coefficients of the radial basis functions with mixed second derivatives. This represents the vertical component of the elastic wave field. In Cartesian coordinate system x The unit vector of direction, These are the coefficients of the vertical second derivative. In Cartesian coordinate system z The unit vector of direction, For the pure transverse wave field vector of the discrete nodes, The transverse wave velocity corresponding to the final discrete node position. i and j All are index numbers. n This represents the total number of nodes participating in the differential calculation.

9. An electronic device, characterized in that, include: One or more processors; A storage device on which one or more programs are stored; When the one or more programs are executed by the one or more processors, the one or more processors cause the one or more processors to implement the method as described in any one of claims 1 to 4.

10. A computer storage medium, characterized in that, It stores a computer program thereon, wherein the computer program, when executed by a processor, implements the method as described in any one of claims 1 to 4.

Images