Novel multi-component conversion multiple wave imaging method for deep-sea rugged seabed structure

Abstract

The invention discloses a novel multi-component conversion multiple imaging method for a deep-sea rugged seabed structure, and relates to the technical field of geophysical exploration for petroleum, and the method comprises the steps: inputting a longitudinal and transverse wave velocity field, a seismic source wavelet, observation system parameters, and a rugged seabed elevation; constructing a non-uniform curved grid of the acoustic-viscoelasticity model; calculating an acoustic-viscoelastic primary wave field continuation operator of vector wave separation under the curved coordinate system; calculating an acoustic-viscoelastic primary wave adjoint wave field of vector wave separation under the curved coordinate system based on an adjoint state theory; calculating an acoustic-viscoelastic n-order multiple wave field continuation operator of vector wave separation under the curved coordinate system; and calculating an acoustic-viscoelastic n-order multiple accompanying wave field of vector wave separation under the curved coordinate system, and generating a multiple imaging result by applying an elastic multiple imaging condition of vector wave separation. According to the invention, full-path compensation and longitudinal and transverse wave vector imaging of multiple waves can be realized, the imaging range is expanded, and the information amount of imaging is increased.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of petroleum geophysical exploration, and in particular to a new method for multi-component converted multiple wave imaging for deep-sea rugged seafloor structure. BACKGROUND

[0002] Seismic data contains primary reflection and multiple waves, and traditional seismic imaging only uses primary reflection, and needs to remove multiple waves in preprocessing to reduce interference. However, multiple waves have long propagation path, wide coverage, carry rich small-angle information and deep structure details, and can provide additional illumination for the underground and expand the imaging range. The existing method realizes imaging by modifying reverse time migration, combining SRME to predict surface multiple waves and inverse scattering series to predict interlayer multiple waves, but faces problems of serious crosstalk noise, low calculation efficiency and poor adaptability to complex geology. Some researches combine least squares reverse time migration with the idea of hierarchical division to preliminarily improve the multiple wave imaging accuracy, but the noise suppression and efficiency still need to be improved.

[0003] Therefore, there is an urgent need for an efficient and accurate multiple wave imaging method to fully utilize multiple wave information, suppress crosstalk noise and improve the imaging quality in complex geological environment. SUMMARY

[0004] To solve the above problems, the present application discloses a new method for multi-component converted multiple wave imaging for deep-sea rugged seafloor structure, which fully utilizes the underground illumination capability of multiple waves and effectively suppresses crosstalk noise to improve the resolution and coverage range of seismic imaging, and is suitable for exploration scenes of complex geological structure.

[0005] To achieve the above purpose, the present application adopts the following technical scheme:

[0006] A new method for multi-component converted multiple wave imaging for deep-sea rugged seafloor structure, comprising the following steps:

[0007] s1. inputting P-S wave velocity field, source wavelet, observation system parameters and rugged seafloor elevation;

[0008] s2. constructing a non-uniform curved grid of acoustic-viscoelastic model;

[0009] s3. calculating acoustic-viscoelastic primary wave wavefield continuation operator of vector wave separation in curved coordinate system;

[0010] s4. calculating acoustic-viscoelastic accompanying wave field of vector wave separation in curved coordinate system based on accompanying state theory;

[0011] s5. calculating acoustic-viscoelastic n-order multiple wave wavefield continuation operator of vector wave separation in curved coordinate system;

[0012] s6. Calculate the vector wavefield of the acoustic-viscoelastic n-order multiple wavefield in the curvilinear coordinate system, and generate the multiple wave imaging result by using the elastic multiple wave imaging condition of the vector wavefield.

[0013] Optionally, in step s2, an acoustic-viscoelastic heterogeneous model grid suitable for the rugged seabed environment is constructed. Seismic waves propagate in the form of acoustic waves in seawater medium in the marine environment, and the first-order acoustic wave equation is used to describe the propagation as shown in the following formula:

[0014] (1);

[0015] wherein, and are the horizontal component and the vertical component of the acoustic wave velocity, P is the acoustic pressure field, f is the source term, x (x, z) is the spatial coordinate, t is the time, λ and μ are the Lame constants, and ρ is the density;

[0016] The relationship between the Lame constant and the velocity is as follows:

[0017] (2);

[0018] wherein, and are the wave velocities of P wave and S wave;

[0019] When the seabed interface is a rugged structure, the spatial coordinates x and z are difficult to accurately obtain. Therefore, the rugged seabed structure model is divided into a non-uniform curvilinear grid, and is mapped into a uniform rectangular grid in the curvilinear coordinate system. After the coordinate transformation, the rugged seabed structure is mapped into a horizontal structure. By using the chain rule, the first-order acoustic wave equation in the curvilinear coordinate system is as follows:

[0020] (3);

[0021] The formula (3) is simplified as:

[0022] (4);

[0023] wherein, is the wavefield continuation operator of the first-order acoustic wave equation in the curvilinear coordinate system, is the acoustic wave field in the fluid medium, P A is the acoustic pressure wave field, and are the horizontal component and the vertical component of the particle velocity, is the source function, and the superscript T is the transpose of the matrix.

[0024] Generally, the deep-sea environment has an irregular seafloor interface, and using finite differences in the traditional Cartesian coordinates (x, z) to image the seafloor structure brings great challenges. To overcome the problem of irregular seafloor, the acoustic-viscoelastic model grid is converted into a curved grid, and coordinate transformation is applied to map the model and equations into a new curved coordinate system (ξ, η), and the equation coefficient matrix , , is calculated by the following formula:

[0025] (5);

[0026] wherein, , , and are the x component and z component of the new curved coordinate system, which are obtained by the following formula:

[0027] (6);

[0028] wherein, (7);

[0029] When the acoustic wave in the fluid layer hits the fluid-solid interface, it will be converted into an elastic wave in the solid layer. Specifically, the pressure field P is converted into the normal stress field and the shear stress field, and various waves are generated in the process, including P waves and converted PS waves. More complexly, the seafloor sediments exhibit viscoelastic behavior, having the characteristics of elastic solids and viscous fluids. To solve the problem of seismic wave propagation in viscoelastic media, the first-order velocity-stress wave equation in viscoelastic isotropic media based on the general standard linear solid (GSLS) model is adopted:

[0030] (8);

[0031] wherein, and are the displacements in the x component and the z component, and are the normal stresses in the x component and the z component, is the shear stress, and are the strain and stress relaxation time of the lth standard linear solid (SLS), v is 1 or 2, depending on the bulk modulus and shear modulus, is 1 or 2, used to distinguish the compression mode and the shear mode; is the number of SLS, v is 1 or 2; , and are memory variables;

[0032] (9);

[0033] and is the non-relaxed Lame constant, and the formula is:

[0034] (10);

[0035] wherein, v is 1 or 2;

[0036] When the seabed interface is a relief interface, the first-order velocity-stress equation of viscoelastic isotropic medium based on the GSLS model in the curvilinear coordinate system is obtained by coordinate transformation:

[0037] (11);

[0038] wherein, and are the relationship coefficients of the strain rate and the stress and the memory variable;

[0039] At the fluid-solid interface, the acoustic wave is converted into the elastic wave, and vice versa; in order to ensure the continuity between the acoustic wave and the elastic wave equation, certain boundary conditions need to be applied, according to:

[0040] (12);

[0041] The two formulas are added:

[0042] (13);

[0043] The boundary condition of the fluid-solid interface is obtained by substituting formula (11) into formula (13):

[0044] (14);

[0045] Formula (14) is valid for both sides of the interface; however, in formula (14), when the fluid-solid interface is a rough terrain, it is difficult to calculate α and β; in addition, due to the stepwise rectangular grid discretization of FDM, there will be artificial false reflections in the conversion process between the acoustic wave and the elastic wave at the fluid-solid interface; when the fluid-solid interface is in a horizontal state, the normal stress and shear stress on the interface disappear;

[0046] The boundary condition of the fluid-solid interface is simplified as:

[0047] (15).

[0048] Alternatively, in step s3, the one-way wave field continuation operator of the acoustic-viscoelastic coupled medium based on vector wave field separation in the curvilinear coordinate system is as follows:

[0049] (16);

[0050] wherein, is the P-wave wavefield extrapolation operator based on wavefield separation of the first-order elastic wave equation in the curvilinear coordinate system, is the S-wave wavefield extrapolation operator based on wavefield separation of the first-order elastic wave equation in the curvilinear coordinate system, is the primary wave wavefield of acoustic wave in fluid medium, , is the primary wave acoustic wavefield in P-wave equation and S-wave equation;

[0051] The acoustic-viscoelastic model grid is converted into a curvilinear grid, and a coordinate transformation is applied to map the model and equation into a new curvilinear coordinate system (ξ, η), and the acoustic wave equation coefficient matrix , , is calculated by the following formula:

[0052] (17);

[0053] wherein, , , and is obtained by the following formula:

[0054] (18);

[0055] wherein, (19);

[0056] In the viscoelastic P-wave equation, , and are the P-wave primary wave particle velocities of x component and z component, and are the P-wave primary wave normal stresses of x component and z component, , and are the mixed primary wave particle velocities of x and z components, and the viscoelastic P-wave wave equation coefficient matrix , , is constructed by the following formula, respectively:

[0057] (20);

[0058] wherein, ;

[0059] is:

[0060] (21);

[0061] wherein, , in which, and are the relaxed Lame constants, and are the number of generalized standard linear solids;

[0062] In the viscoelastic shear wave equation, , and are the shear wave primary wave particle velocity, and are the shear wave primary wave normal stress, is the shear wave primary wave shear stress, , the viscoelastic shear wave wave equation coefficient matrix , , is:

[0063] (22);

[0064] wherein, is:

[0065] (23);

[0066] wherein, is the memory variable;

[0067] The relationship between the relaxed and non-relaxed Lame constants is:

[0068] (24);

[0069] Mixed velocity wavefield ( , ), P-wave field ( , ) and S-wave field ( , ) satisfy:

[0070] (25);

[0071] , , and are updated by:

[0072] (26);

[0073] wherein, , , , the coefficients J1, J2 are given by the following formula:

[0074] (27);

[0075] wherein, .

[0076] Optionally, in step s4, the adjoint wavefield of the acoustic-elastic coupling one-way wave receiver point back-propagation based on vector wavefield separation of the adjoint state theory in the curvilinear coordinate system is calculated by the following formula:

[0077] (28);

[0078] wherein, is the adjoint matrix of matrix , is the adjoint wavefield of , , and are the acoustic wave, viscoelastic longitudinal wave and viscoelastic transverse wave adjoint operators in the curvilinear coordinate system respectively, , is the η coordinate value at the seabed interface, is the maximum η coordinate value, is the maximum acoustic wave calculation depth in the wavefield back-propagation process, and M is the finite difference precision, and are calculated by the following formula respectively:

[0079] (29);

[0080] wherein, is the synthetic acoustic pressure record, is the observed acoustic pressure record, and are the observed and synthetic elastic displacement fields respectively, and the subscripts x and z are the horizontal and vertical component displacements respectively, and are the mixed records after P-wave and S-wave wavefield separation, and are obtained by the Born approximation linear forward method.

[0081] Optionally, in step s5, based on the Born approximation theory, the first arrival wave synthetic data of the coupled acoustic wave-decoupled elastic wave is represented as follows:

[0082] (30);

[0083] where, , , , and denote the perturbation wavefield of the primary wave , , , , and x r denotes the position of the receiver, denotes the minimum elastic calculation depth in the wavefield inversion process, and denote the synthetic P-wave and S-wave records of the first arrival wave, , and denote the reconstructed acoustic, P-wave, and S-wave source terms of the first arrival wave, which are calculated by the following formula:

[0084] (31);

[0085] where, and denote the first arrival wave imaging results of the P-wave and S-wave components, respectively;

[0086] Similarly, the synthetic data of the first-order, second-order, …, n-order multiple waves are calculated in the following manner:

[0087] (32);

[0088] where, , , , and denote the perturbation wavefield of the first-order multiple wave, , , , and denote the perturbation wavefield of the n-order multiple wave, and denote the synthetic P-wave and S-wave records of the first-order multiple wave, and denote the synthetic P-wave and S-wave records of the n-order multiple wave; the reconstructed acoustic, P-wave, and S-wave source terms of the n-order multiple wave , and are calculated by the following formula:

[0089] (33);

[0090] where, and respectively represent the imaging results of the n-th order wave of the longitudinal wave and the transverse wave component.

[0091] Optionally, in step s6, the elastic multi-component n-th order multiple wave imaging condition of the vector wave separation is:

[0092] (34);

[0093] wherein, and are the imaging results of the n-th order multiple wave of the longitudinal wave and the transverse wave, , , , are the n-th order multiple wave fields of the x-component longitudinal wave, the x-component transverse wave, the z-component longitudinal wave and the z-component transverse wave respectively, , , , are the n-th order multiple wave accompanying wave fields of the x-component longitudinal wave, the x-component transverse wave, the z-component longitudinal wave and the z-component transverse wave respectively, and are obtained by the following formula:

[0094] (35);

[0095] wherein, , and are the accompanying wave fields of the n-th order multiple wave fields , and , is the maximum acoustic wave calculation depth in the process of wave field reverse propagation, and p n represents the observed acoustic pressure record of the first arrival wave (n=0) and the n-th order multiple wave (n>0), is the observed elastic seismic record representing the first arrival wave (n=0) and the n-th order multiple wave (n>0).

[0096] The beneficial effects of the present application are that the present application innovatively constructs a full-path propagation operator of an elastic surface layer multiple wave field, and through a vector wave field separation method, proposes a sound-elastic coupling surface layer multiple wave imaging technology for a fluctuating seabed. The method breaks through the limitations of traditional elastic primary wave imaging, and compared with single-component multiple wave imaging, it can provide more abundant seismic imaging information. Specifically, the technology fully utilizes the multiple scattering and conversion effects of the complex structure of the seabed on the wave field, so that different types of wave field information can be effectively extracted and integrated, thereby improving the clarity and accuracy of imaging. The fluctuating seabed sound-elastic coupling multiple wave imaging method provides more abundant seismic imaging information than conventional elastic primary wave imaging and single-component multiple wave imaging. BRIEF DESCRIPTION OF DRAWINGS

[0097] Figure 1A flowchart of a new method for multi-component converted multiple wave imaging facing deep-sea rugged seafloor structure according to the present application is shown in the figure;

[0098] Figure 2 A Sigsbee2B model graph of a fluctuating seafloor acoustic medium according to an embodiment of the present application is shown in the figure;

[0099] Figure 3 A Sigsbee2B reflection coefficient model graph of an acoustic medium according to an embodiment of the present application is shown in the figure;

[0100] Figure 4 A mesh partitioning schematic diagram according to an embodiment of the present application is shown in the figure;

[0101] Figure 5 An elastic primary wave shot record of wavefield separation according to an embodiment of the present application is shown in the figure;

[0102] Figure 6 An elastic first-order multiple wave shot record of wavefield separation according to an embodiment of the present application is shown in the figure;

[0103] Figure 7 An elastic second-order multiple wave shot record of wavefield separation according to an embodiment of the present application is shown in the figure;

[0104] Figure 8 A primary wave imaging graph according to an embodiment of the present application is shown in the figure;

[0105] Figure 9 Multiple wave imaging according to an embodiment of the present application is shown in the figure. DETAILED DESCRIPTION

[0106] In order to make the objects, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the accompanying drawings in the embodiments of the present application. Obviously, the described embodiments are only some but not all of the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative work fall within the scope of protection of the present application. Therefore, the following detailed description of the embodiments of the present application provided in the accompanying drawings is not intended to limit the scope of the claimed present application, but only for selected embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative work fall within the scope of protection of the present application.

[0107] A new method for multi-component converted multiple wave imaging facing deep-sea rugged seafloor structure, as shown in the figure, comprises the following steps: Figure 1

[0108] ​s1. input P-S wave velocity field, source wavelet, observation system parameters, rugged seafloor elevation.

[0109] s2. construct non-uniform curved grid of acoustic-viscoelastic model; construct acoustic-viscoelastic non-uniform model grid suitable for rugged seafloor environment, seismic wave first propagates in the form of acoustic wave in seawater medium in marine environment, and is described by the following first-order acoustic wave equation:

[0110] (1);

[0111] wherein, and are horizontal component and vertical component acoustic wave velocity respectively, P is acoustic pressure field, f is source term, x (x, z) is spatial coordinate, t is time, λ and μ are Lame constants, and ρ is density;

[0112] The relationship between Lame constant and velocity is:

[0113] (2);

[0114] wherein, and are P-wave and S-wave velocity respectively;

[0115] When the seafloor interface is a severe rugged structure, the spatial coordinates x and z are difficult to accurately obtain, therefore, the rugged seafloor structure model is divided into non-uniform curved grid, and is mapped into uniform rectangular grid in curved coordinate system, after coordinate transformation, the rugged seafloor structure is mapped into horizontal structure; by applying chain rule, the first-order acoustic wave equation in curved coordinate system is:

[0116] (3);

[0117] Equation (3) is simplified as:

[0118] (4);

[0119] wherein, is wave field continuation operator of the first-order acoustic wave equation in curved coordinate system, is acoustic wave field in fluid medium, P A is acoustic pressure wave field, and are horizontal component and vertical component particle velocity respectively, is source function, and superscript T is transpose of matrix;

[0120] Generally, the deep-sea environment has an irregular seafloor interface, and using finite differences in the traditional Cartesian coordinates (x, z) to image the seafloor structure poses a great challenge. To overcome the problem of the irregular seafloor, the acoustic-viscoelastic model grid is converted into a curvilinear grid, and a coordinate transformation is applied to map the model and equations into a new curvilinear coordinate system (ξ, η). The coefficient matrix of the equation is , , calculated by the following formula:

[0121] (5);

[0122] wherein , , and are the x and z components of the new curvilinear coordinate system, which are obtained by the following formula:

[0123] (6);

[0124] wherein (7);

[0125] When the acoustic wave in the fluid layer hits the fluid-solid interface, it will be converted into an elastic wave in the solid layer. Specifically, the pressure field P is converted into a normal stress field and a shear stress field, in the process of which various waves are generated, including P waves and converted PS waves. More complexly, the seafloor sediments exhibit viscoelastic behavior, having the characteristics of elastic solids and viscous fluids. To solve the problem of the propagation of seismic waves in viscoelastic media, a first-order velocity-stress wave equation in viscoelastic isotropic media based on the general standard linear solid (GSLS) model is adopted:

[0126] (8);

[0127] wherein and are the displacements in the x and z components, and are the normal stresses in the x and z components, is the shear stress, and are the strain and stress relaxation time of the lth standard linear solid (SLS), v is 1 or 2, depending on the bulk modulus and shear modulus, is 1 or 2, used to distinguish the compression mode and the shear mode; is the number of SLS, v is 1 or 2; , and are memory variables;

[0128] (9);

[0129] and is the unrelaxed Lame constant, and the formula is:

[0130] (10);

[0131] where, v is 1 or 2;

[0132] When the seabed interface is a relief interface, the first-order velocity-stress equation of viscoelastic isotropic medium based on the GSLS model in the curvilinear coordinate system is obtained by coordinate transformation:

[0133] (11);

[0134] where, and are the relationship coefficients of strain rate and stress, and memory variable;

[0135] At the fluid-solid interface, the sound wave is converted into an elastic wave, and vice versa; in order to ensure the continuity between the sound wave and the elastic wave equation, certain boundary conditions need to be applied, according to:

[0136] (12);

[0137] Add the two equations:

[0138] (13);

[0139] Substitute equation (11) into equation (13) to obtain the boundary condition of the fluid-solid interface:

[0140] (14);

[0141] Equation (14) is valid for both sides of the interface; however, in equation (14), when the fluid-solid interface is a rough terrain, it is difficult to calculate α and β; in addition, due to the stepwise rectangular grid discretization of FDM, there will be artificial false reflections in the conversion process between the sound wave and the elastic wave at the fluid-solid interface; when the fluid-solid interface is in a horizontal state, the normal stress and shear stress on the interface disappear;

[0142] The boundary condition of the fluid-solid interface is simplified as:

[0143] (15).

[0144] The present application discretizes fluid layers (acoustic medium) and solid layers (viscoelastic medium) into curved meshes, and converts them into corresponding rectangular meshes. Based on this strategy, irregular fluid-solid interfaces can be converted into horizontal interfaces. The acoustic wave equation, the elastic wave equation and the fluid-solid boundary conditions in the physical domain (x,z) are converted into conditions in the computational domain (ξ,η).

[0145] s3. One-way wavefield extrapolation operator for acoustic-viscoelastic coupling media based on vector wavefield separation in curvilinear coordinate system As follows:

[0146] (16);

[0147] wherein, is a P-wave wavefield extrapolation operator for the first-order elastic wave equation based on wavefield separation in the curvilinear coordinate system, is an S-wave wavefield extrapolation operator for the first-order elastic wave equation based on wavefield separation in the curvilinear coordinate system, is an acoustic one-way wavefield in a fluid medium, , is a one-way acoustic wavefield in the P-wave equation and the S-wave equation;

[0148] The acoustic-viscoelastic model mesh is converted into a curved mesh, and a coordinate transformation is applied to map the model and the equation into a new curved coordinate system (ξ,η). The acoustic wave equation coefficient matrix , , is calculated by the following formula:

[0149] (17);

[0150] wherein, , , and are obtained by the following formula:

[0151] (18);

[0152] wherein, (19);

[0153] In the viscoelastic P-wave equation, , and are the P-wave one-way particle velocities in the x component and the z component, and are the P-wave one-way normal stresses in the x component and the z component, , and mixed P-wave particle velocities in x and z components, respectively, coefficient matrix of viscoelastic longitudinal wave equation , , are constructed by the following equations, respectively:

[0154] (20);

[0155] wherein, ,

[0156] (21);

[0157] wherein, , in the equation, and are relaxed Lame constants, and are the number of generalized standard linear solids;

[0158] In the viscoelastic transverse wave equation, , and are transverse wave particle velocities, and are transverse wave primary wave normal stresses, is a transverse wave primary wave shear stress, coefficient matrix of viscoelastic transverse wave equation , , are:

[0159] (22);

[0160] wherein, are:

[0161] (23);

[0162] wherein, is a memory variable;

[0163] The relationship between the relaxed and non-relaxed Lame constants is:

[0164] (24);

[0165] Mixed velocity wave field ( , ), P-wave field ( , ) and S-wave field ( , ​) are satisfied:

[0166] (25);

[0167] , , and are updated by

[0168] (26);

[0169] where, , , , the coefficients J1, J2 are given by

[0170] (27);

[0171] where, .

[0172] s4. Calculate the acoustic-viscoelastic primary wave accompanying wave field of vector wave separation in the curvilinear coordinate system based on the accompanying state theory , which is calculated by

[0173] (28);

[0174] where, is the adjoint matrix of matrix , is the adjoint wave field of , , and , and are the acoustic wave, viscoelastic longitudinal wave and viscoelastic transverse wave accompanying operators in the curvilinear coordinate system respectively, , is the η coordinate value at the seabed interface, is the maximum η coordinate value, is the maximum acoustic wave calculation depth in the process of wave field reverse propagation, and M is the finite difference precision, and are calculated by

[0175] (29);

[0176] where, is the synthetic acoustic pressure record, is the observed acoustic pressure record, and are the observed and synthesized elastic displacement fields respectively, and the subscripts x and z are the horizontal and vertical component displacements respectively, and mixed records that are separated by P-wave and S-wave wavefield, and are obtained by the Born approximation based linear forward modeling method.

[0177] s5. Computing the acoustic-viscoelastic n-order multiple wavefield continuation operator of vector wave separation in the curvilinear coordinate system where the first arrival wave synthetic data that couples acoustic wave-decouples elastic wave is represented as follows based on the Born approximation theory:

[0178] (30);

[0179] where, , , , and represent the perturbation wavefield of the first-order wave , , , , and x r represents the position of the receiver, represents the minimum elastic calculation depth in the wavefield inversion process, and represent the synthetic P-wave and S-wave records of the first arrival wave, , and represent the reconstructed acoustic wave, P-wave, and S-wave source terms, which are calculated by the following formula:

[0180] (31);

[0181] where, and represent the first arrival wave imaging results of the P-wave and S-wave components, respectively;

[0182] Similarly, the synthetic data of the first-order, second-order, …, n-order multiple waves are calculated as follows:

[0183] (32);

[0184] where, , , , and represent the perturbation wavefield of the first-order multiple wave, , , , and represent the perturbation wavefield of the n-order multiple wave, and representing the synthetic P and S records of the first-order multiples, and representing the synthetic P and S records of the n-order multiples; the reconstructed acoustic, P and S source terms of the n-order multiples , and are calculated by

[0185] (33).

[0186] where, and represent the imaging results of the n-order wave of P and S components, respectively.

[0187] s6. Calculate the acoustic-viscoelastic n-order multiple accompanying wave field of the vector wave separation in the curvilinear coordinate system, and generate the multiple imaging results by applying the vector wave separation elastic multiple imaging conditions; the vector wave separation elastic multiple n-order multiple imaging conditions are:

[0188] (34).

[0189] where, and are the imaging results of the n-order multiple of P and S waves, , , , are the n-order multiple wave fields of the x-component P wave, the x-component S wave, the z-component P wave, and the z-component S wave, respectively, , , , are the n-order multiple accompanying wave fields of the x-component P wave, the x-component S wave, the z-component P wave, and the z-component S wave, respectively, and are calculated by

[0190] (35).

[0191] where, , and are the accompanying wave fields of the n-order multiple wave fields , and , is the maximum acoustic calculation depth in the wave field reverse propagation process, and p n represents the observed acoustic pressure record of the first arrival and n-order multiple, is the observed elastic seismic record representing the first arrival and n-order multiple.

[0192] Test example

[0193] To verify the effectiveness of the proposed method, the acoustic-elastic multiple wave simulation is performed on the Sigsbee2B model. This model is widely used in the field of marine exploration and contains complex subsurface geological structures, which can be used to test the performance of multiple wave imaging methods in actual complex geology, such as Figure 2 As shown in the figure, the left figure is the P-wave velocity model, and the right figure is the S-wave velocity model. Figure 3 The corresponding reflection coefficient model of the model is shown. The reflection coefficient describes the acoustic-elastic impedance change between different subsurface media and is one of the key parameters for reflection wave simulation. The reflection coefficient model provides the necessary physical background for multiple wave propagation and helps to accurately simulate the multiple wave propagation characteristics of the subsurface interface and its complex structure. In addition, in order to improve the accuracy and computational efficiency of numerical simulation, the body-fitted grid partitioning technique is used to discretize the model. Figure 4 The grid partitioning is fine and closely matches the geometry of the subsurface medium, ensuring high accuracy of the simulation results. When simulating multiple wave propagation, the accuracy of grid partitioning directly affects the propagation process of reflected and multiple waves. Fine grid can accurately capture the small differences between strata, further improving the accuracy of reflected and multiple wave imaging of the model. Figure 5 The wave field separated elastic primary wave shot record obtained by forward simulation is shown, which includes P-wave horizontal component, S-wave horizontal component, P-wave vertical component, and S-wave vertical component. From these images, it can be seen that the propagation path of elastic primary wave in x and z components is clear, and the propagation characteristics of P-wave and S-wave are accurately reproduced.

[0194] Figure 6 The wave field separated elastic first-order multiple wave shot record obtained by forward simulation is shown, which includes P-wave horizontal component, S-wave horizontal component, P-wave vertical component, and S-wave vertical component. It can be seen that the reflection and refraction processes experienced by the first-order multiple wave during propagation. Although the propagation path of multiple wave is more complex than that of primary wave, the P-wave and S-wave components of multiple wave are accurately simulated in all components. These results show that the method not only can effectively simulate the propagation process of primary wave, but also can successfully handle the propagation characteristics of different order multiple waves, providing deep insights into the multiple wave propagation process.

[0195] Figure 7The elastic second-order multiple wave shot record separated by the wave field obtained by forward simulation is in turn: P-wave horizontal component, S-wave horizontal component, P-wave vertical component, and S-wave vertical component, which shows the propagation path of the second-order multiple wave and the complexity of its P-wave and S-wave components. Although the propagation path of the second-order multiple wave is more tortuous, and the reflection and refraction mechanisms are more complex, the forward simulation method can still accurately reproduce the waveform characteristics, especially in the vertical component and the horizontal component, the distribution of P-wave and S-wave is ideally simulated, which further verifies the efficiency and accuracy of the forward simulation method in processing high-order multiple waves, and provides reliable technical support for wave field simulation in complex geological environments.

[0196] Further imaging is performed using the method of the application and compared with the results of conventional imaging methods. Figure 8 The imaging results of the elastic primary wave are shown, the left image is the longitudinal wave imaging, and the right image is the transverse wave imaging; it can be seen that the primary wave imaging shows the basic characteristics of the underground structure in the performance of longitudinal and transverse waves, and can effectively identify the velocity change and shear characteristics of the medium. However, due to the limitations of the primary wave in the propagation process in the complex underground medium, the imaging resolution and coverage range are relatively limited, especially in some more complex geological structures, the imaging results may appear certain shadow or distortion phenomenon. Figure 9 The imaging results of the elastic multiple wave are shown, the left image is the longitudinal wave imaging, and the right image is the transverse wave imaging. Compared with the primary wave imaging, the multiple wave imaging shows more clear and detailed underground structure details, especially in high-resolution imaging. Through the propagation path of the multiple wave, the imaging method can effectively penetrate the complex geological layers underground, make up for the blind area in the primary wave imaging, and significantly improve the spatial resolution of the imaging. In addition, the reflection characteristics of the multiple wave enable it to cover a wider area, especially in the areas that cannot be fully detected by traditional primary wave imaging, successfully achieving more extensive scanning and analysis of the underground structure.

[0197] By comparing the results of the primary wave imaging and the multiple wave imaging, it can be obviously found that the proposed multiple wave imaging method has significant advantages in resolution and imaging range. The multiple wave imaging can clearly reveal the details of the complex underground structure, especially in the areas that the primary wave cannot effectively penetrate, providing more comprehensive geological information; at the same time, the multiple wave imaging can cover a wider area, expanding the application range of the traditional imaging method. Therefore, the imaging method based on the multiple wave has higher resolution and wider imaging range in actual seismic exploration, providing a more accurate and comprehensive solution for exploration in complex geological environments.

[0198] Of course, the above description is not a limitation of the present application, and the present application is not limited to the above examples. Changes, modifications, additions or substitutions made by those skilled in the art within the spirit and scope of the present application should also be included in the protection scope of the present application.

Claims

1. A novel multi-component converted multiple-wave imaging method for deep-sea rugged seabed structures, characterized in that, Includes the following steps: s1. Input the P-wave and S-wave velocity fields, source wavelet, observation system parameters, and rugged seabed elevation; s2. Construct a non-uniform curved mesh for the acoustic-viscoelastic model; s3. Calculate the acoustic-viscoelastic primary wave field extension operator for vector wave separation in curved coordinates; s4. Calculate the adjoint wave field of the acoustic-viscoelastic primary wave in curved coordinates based on the adjoint state theory; s5. Calculate the acoustic-viscoelastic nth-order wavefield extension operator for vector wave separation in curved coordinates; s6. Calculate the accompanying wave field of the acoustic-viscoelastic nth-order multiple wave separated by vector wave in curved coordinate system, and generate multiple wave imaging results by applying the elastic multiple wave imaging conditions of vector wave separation.

2. The novel multi-component converted multiple-wave imaging method for deep-sea rugged seabed structures as described in claim 1, characterized in that, In step s2, an acoustic-viscoelastic non-uniform model mesh suitable for rugged seabed environments is constructed. Seismic waves first propagate in the seawater medium as sound waves in the marine environment, described by the following first-order acoustic wave equation: （1）； in, and Let x be the sound wave velocity of the horizontal component and the vertical component, respectively; P be the sound pressure field; f be the source term; x(x,z) be the spatial coordinates; t be the time; λ and μ be the Lamé constants; and ρ be the density. The relationship between Lamé constant and velocity is as follows: （2）； in, and These are the wave velocities of P-waves and S-waves, respectively. The rugged seafloor structure model is divided into a non-uniform curved mesh and mapped to a uniform rectangular mesh in a curved coordinate system. After coordinate transformation, the rugged seafloor structure is mapped to a horizontal structure. Applying the chain rule, the first-order acoustic wave equation in the curved coordinate system is: （3）； Equation (3) can be simplified as follows: （4）； in, For the wavefield extension operator of the first-order acoustic wave equation in curved coordinates, For the acoustic wave field in the fluid medium, P A For sound pressure wave field, and These are the particle velocities of the horizontal and vertical components, respectively. Here is the source function, and the superscript T is the transpose of the matrix; The acoustic-viscoelastic model is meshed into a curved mesh, and a coordinate transformation is applied to map the model and equations to a new curved coordinate system (ξ, η). The equation coefficient matrix... , , It is calculated by the following formula: （5）； in, , , and The x and z components of the new curvilinear coordinate system are obtained from the following equation: （6）； in, (7); To address the propagation problem of seismic waves in viscoelastic media, the first-order velocity-stress wave equation for viscoelastic isotropic media based on the general standard linear solid GSLS model is adopted: （8）； in, and Let x be the displacements in the x and z components. and Let x and z be the normal stresses, respectively, for the x and z components. For shear stress, and ν represents the strain and stress relaxation time of the l-th standard linear solid SLS, respectively, where v is 1 or 2, depending on the bulk modulus and shear modulus. It can be either 1 or 2, used to distinguish between compression mode and shearing mode; The SLS number, where v is 1 or 2; , and It is a memory variable; （9）； and The non-relaxed Lamé constant is calculated using the following formula: （10）； in, v is 1 or 2; When the seabed interface is undulating, a coordinate transformation is performed to obtain the first-order velocity-stress equation for the viscoelastic isotropic medium based on the GSLS model in curved coordinates: （11）； in, and The coefficients representing the relationship between strain rate and stress, and memory variables; At the fluid-solid interface, sound waves are converted into elastic waves, and vice versa; to ensure the continuity between the equations for sound waves and elastic waves, according to: （12）； Add these two expressions together: （13）； Substituting equation (11) into equation (13), we obtain the boundary conditions for the fluid-solid interface: （14）； Equation (14) is valid on both sides of the interface; The boundary conditions at the fluid-solid interface are simplified as follows: （15）。 3. A novel multi-component converted multiple-wave imaging method for deep-sea rugged seabed structures as described in claim 1, characterized in that, In step s3, the primary wave field extension operator of the acoustic-viscoelastic coupling medium based on vector wave field separation in curved coordinates is used. as follows: （16）； in, For the longitudinal wave field extension operator of the first-order elastic wave equation based on wave field separation in curved coordinate system, This is the transverse wavefield extension operator for the first-order elastic wave equation based on wavefield separation in curved coordinates. The first wave field of sound in a fluid medium. , The primary wave field of the P-wave equation and the S-wave equation; The acoustic-viscoelastic model is meshed into a curved mesh, and a coordinate transformation is applied to map the model and equations to a new curved coordinate system (ξ, η). The coefficient matrix of the acoustic wave equation is... , , It is calculated by the following formula: （17）； in, , , and It can be obtained from the following formula: （18）； in, (19); In the viscoelastic longitudinal wave equation , and The values ​​represent the particle velocities of the primary longitudinal wave in the x and z components, respectively. and The primary wave normal stresses of the longitudinal wave, representing the x-component and the z-component, are respectively. , and The mixed primary wave particle velocities of the x and z components, respectively, and the coefficient matrix of the viscoelastic longitudinal wave equation. , , They are constructed from the following formulas respectively: （20）； in, ; for: （21）； in, In the formula, and Let Lamé be the relaxed constant. and The number of generalized standard linear solids; In the viscoelastic transverse wave equation , and The velocity of the primary wave particles in the transverse wave. and For the normal stress of the first wave of the transverse wave, For the primary wave shear stress of the transverse wave, viscoelastic transverse wave equation coefficient matrix , , for: （22）； in, for: （23）； in, For memory variables; The relationship between the relaxed and unrelaxed Lamé constants is as follows: （24）； Mixed velocity wave field ( , P-wave field ( , ) and S-wave field ( , )satisfy: （25）； , , and The following formula is updated to obtain: （26）； in, , , The coefficients J1 and J2 are given by the following formula: （27）； in, .

4. A novel multi-component converted multiple-wave imaging method for deep-sea rugged seabed structures as described in claim 1, characterized in that, In step s4, the accompanying wave field of the acoustic-elastic coupling primary wave detector backpropagation based on the adjoint state theory in the curved coordinate system is obtained by vector wave field separation. It is calculated by the following formula: （28）； in, For matrix The adjoint matrix, for The accompanying wave field, , and These are the adjoint operators for sound waves, viscoelastic longitudinal waves, and viscoelastic transverse waves in curvilinear coordinate systems, respectively. , Let η be the coordinate value at the seabed interface. For the maximum η coordinate value, Let M be the maximum acoustic wave depth calculated during the back-propagation process in the wave field, and M be the finite difference accuracy. and Calculated using the following formulas respectively: （29）； in, For synthesized sound pressure recording, For the observed sound pressure records, and These are the observed and synthesized elastic displacement fields, respectively, with subscripts x and z representing the horizontal and vertical displacement components, respectively. and This is a mixed record after P-wave and S-wave field separation. and It was obtained using the Born approximate linear forward modeling method.

5. A novel multi-component converted multiple-wave imaging method for deep-sea rugged seabed structures as described in claim 4, characterized in that, In step s5, based on the Born approximation theory, the first arrival synthesized data of the coupled acoustic wave and decoupled elastic wave are represented as follows: （30）； in, , , , and Indicates a first wave , , , ,and The perturbation wave field, x r Indicates the location of the detector point. This represents the minimum elastic calculation depth during wavefield inversion. and This represents the composite P-wave and S-wave record of the first arrival wave. , and The reconstructed acoustic, P-wave, and S-wave source terms are calculated using the following formula: （31）； in, and These represent the first arrival imaging results for the P-wave and S-wave components, respectively. Similarly, the composite data of first-order, second-order...n-order multiple waves are calculated as follows: （32）； in, , , , and This represents the perturbation wave field of a first-order multiple wave. , , , and This represents the perturbation wave field of an nth-order multiple wave. and This represents the composite P-wave and S-wave record of a first-order multiple wave. and Represents the synthetic P-wave and S-wave records of an nth-order multiple; the source terms of the acoustic, P-wave, and S-wave of the reconstructed nth-order multiple. , and It is calculated by the following formula: （33）； in, and These represent the imaging results of the nth-order waves of the longitudinal and transverse wave components, respectively.

6. A novel multi-component converted multiple-wave imaging method for deep-sea rugged seabed structures as described in claim 1, characterized in that, In step s6, the conditions for elastic multi-component n-order multiple imaging of vector wave separation are: （34）； in, and The imaging results for the longitudinal and transverse waves of an nth-order multiple wave are shown. , , , These are the nth-order wavefields of the x-component longitudinal wave, x-component transverse wave, z-component longitudinal wave, and z-component transverse wave, respectively. , , , The accompanying wave fields of the nth-order multiple waves, representing x-component P-waves, x-component S-waves, z-component P-waves, and z-component S-waves respectively, are obtained by the following equation: （35）； in, , and For nth-order multiple wave fields , and The accompanying wave field, The maximum acoustic wave computational depth during the back-propagation of the wave field, p n This represents the observed sound pressure records of the first arrival wave and the nth-order multiple wave. This represents the observed elastic seismic records of the first arrival wave and the nth-order multiple waves.

Images