An Application of 3D Cross-Well Elastic Reverse Time Migration Imaging Based on the Multi-Wave and Multi-Component Technique in Coastal Engineering Exploration

: Precise surveys are indispensable in coastal engineering projects. The extensive presence of sand in the coastal area leads to significant attenuation of seismic waves within unsaturated loose sediments. As a result, it becomes challenging for seismic waves to penetrate the weathered zone and reach the desired depth with significant amount of energy. In this study, the application of three-dimensional (3D) cross-well elastic reverse time migration (RTM) imaging based on multi-wave and multi-component techniques in coastal engineering exploration is explored. Accurate decomposition of vector compressional (P) and shear (S) waves is achieved through two wavefield decoupling algorithms without any amplitude and phase distortion. Additionally, compressional wave pressure components are obtained, which facilitates subsequent independent imaging. This study discusses and analyzes the imaging results of four imaging strategies under cross-correlation imaging conditions in RTM imaging. The analysis leads to the conclusion that scalarizing vector wavefields imaging yields superior imaging of P-and S-waves. Furthermore, the imaging results obtained through this approach are of great physical significance. In order to validate the efficacy of this method in 3D geological structure imaging in coastal areas, RTM imaging experiments were performed on two representative models. The results indicate that the proposed 3D elastic wave imaging method effectively generates accurate 3D cross-well imaging of P-and S-waves. This method utilizes the multi-wave and multi-component elastic wave RTM imaging technique to effectively leverage the Earth’s elastic medium without increasing costs. It provides valuable information about the distribution of subsurface rock layers, interfaces, and other structures in coastal engineering projects. Importantly, this can be achieved without resorting to extensive excavation or drilling operations. This method addresses the limitations of current cross-well imaging techniques, thereby providing abundant and accurate geological and geophysical information for the analysis and interpretation of 3D geological structures in coastal engineering projects. It has important theoretical and practical significance in real-world production, as well as for the study of geological structures in coastal engineering.


Introduction
During coastal engineering projects, extensive surveys and investigations are crucial for gathering information on geological structures, sediment distribution, and soil properties in coastal areas.These investigations involve analyses by drilling, sampling, and geophysical exploration to evaluate the geological formations, soil types, lithology, and groundwater levels [1][2][3].Among the numerous geophysical exploration methods, seismic exploration is widely employed due to its ability to provide superior subsurface resolution [4].However, the extensive presence of sand in the coastal area leads to a significant attenuation of seismic waves within unsaturated loose sediments.Surface seismic methods are highly susceptible to variations in near-surface conditions and encountering difficulties in transmitting a significant amount of energy through the weathered zone to reach the targeted depths [5][6][7].Compared to surface seismic reflection or refraction surveys, cross-well seismic acquisition methods emerge as an optimal choice for effectively targeting specific coastal areas [8,9].
Cross-well seismic can be applied to various aspects of oil and gas exploration, including detailed imaging of structures and precise characterization of reservoirs [10][11][12].Furthermore, cross-well seismic has also been rapidly developed in the field of engineering [13,14] and has been applied in diverse areas, such as geological engineering [15][16][17], hydrogeological surveys [18][19][20], and quality inspections in civil engineering projects [21][22][23].Currently, most cross-well seismic studies focus on two-dimensional (2D) tomographic imaging between adjacent wells [24,25].It requires data collection and travel-time inversion between two wells to acquire subsurface structural profiles.Nevertheless, this imaging approach has certain limitations.For example, in 2D computed tomography (CT) imaging, only adjacent wells can be imaged, and information from multiple wells distributed in 3D space cannot be effectively utilized [26].As a result, the utilization rate of spatial well locations remain relatively low.Furthermore, the 2D cross-well CT imaging technique can only provide information on the geological structure of the profile between the two wells, failing to capture the lateral structure information of the profile [27].As a result, the imaging provides poor lateral continuity, hindering the evaluation of 3D geological structures.
Additionally, there are issues associated with imaging steep-dip interfaces (interfaces with dip angles exceeding 45 • ) [28].In recent years, the rapidly developing RTM method follows the full-wave wave equation during wavefield extrapolation and is not limited by angles [29][30][31].Among various pre-stack migration methods, RTM is considered the most accurate imaging technique.However, current studies on cross-well RTM imaging are limited to 2D space or 3D P-waves RTMs [32,33].With the advancement of seismic exploration technology, the challenges have become increasingly complex.Conventional compressional wave exploration can no longer satisfy the needs of engineering projects [34].Improvements need to be made in two key aspects: First, it is necessary to image areas with weak compressional wave reflection energy or no compressional wave reflection signals to provide reliable structural imaging information.Second, it is necessary to provide imaging information for elastic multi-wave analysis, thus establishing a reliable foundation for subsequent elastic wave amplitude analysis and rock physics parameter inversion.Conventional seismic exploration techniques commonly employ compressional wave sources.To address the limitations of compressional wave reflection energy, the shear waves generated through compressional wave conversion need to be imaged [35,36].This approach aims to fully utilize the Earth's elastic media to construct comprehensive imaging information of subsurface structures without increasing costs.Simultaneously, the imaging results must include multi-wave imaging results for subsequent amplitude variation with angle (AVA) analysis and the inversion of rock physics parameters [37,38].Compared with the conventional compressional wave exploration techniques, the theoretical assumptions of multi-wave and multi-component seismic exploration technology are more consistent with the actual characteristics of subsurface media, offering significant advantages in addressing practical challenges.
To effectively address the above issues, it is essential to fully utilize the distinctive features of multi-wave and multi-component seismic data acquired between wells.Therefore, this paper investigates the application of 3D cross-well elastic wavefield RTM imaging in coastal engineering exploration based on multi-wave and multi-component data.By applying multi-wave and multi-component elastic wave RTM imaging techniques, precise 3D geological structures between wells can be obtained, laying the foundation for subsequent lithological analysis and inversion.This approach addresses the limitations of current cross-well imaging techniques, thus providing rich geological and geophysical information for analyzing and interpreting 3D geological structures in coastal regions.It has important theoretical and practical significance in actual production, as well as for the investigation of geological structures in coastal engineering.

Methodology
To study the application of 3D cross-well seismic exploration in coastal engineering, it is essential to first conduct research on the numerical simulation of elastic waves.Numerical simulations encompass the simulation of the propagation characteristics of elastic waves in various media through mathematical methods based on the known physical parameters of the detection model, combined with elastic wave propagation theory [39].Numerical simulations can be divided into three categories: wave equation numerical methods, integral methods, and ray tracing methods.Among these, wave equation numerical methods primarily include the finite difference method (FDM) [40], pseudo spectral method (PSM) [41], finite element method (FEM) [42], spectral element method (SEM) [43], and boundary element method (BEM) [44].FDM is extensively employed in numerical simulations of elastic waves due to its advantageous features, including fast computation speed and high accuracy.
In this section, we begin by discussing the 3D elastic wave equations of motion and grid discretization in cylindrical coordinates.Then, we utilize two different methods for elastic wave decomposition.Following this, we introduce the principles of RTM imaging and provide detailed explanations of four different cross-correlation imaging conditions.Lastly, we propose the Poynting vector and Laplace filtering to attenuate RTM imaging artifacts.

3D Elastic Wave Equations of Motion and Grid Discretization
The basic laws governing the propagation of elastic waves are described by the wave equation of elastic waves.Based on the equations of motion, Cauchy equations, and physical equations in a three-dimensional Cartesian coordinate system, the first-order velocity-stress equation in an isotropic medium in cylindrical coordinates can be expressed as follows [45]: where v = v x , v y , v z T denotes the velocity; τ = τ xx , τ yy , τ zz , τ xy , τ yz , τ xz T denotes the stress vectors; f = f x , f y , f z T represents the point force source; g = g xx , g yy , g zz , g xy , , g xz , g yz T represents the coupling; ρ is the density; λ and µ refer to the Lamé constants.

Elastic Wave Decomposition
For elastic wave RTM, the key to obtaining multi-wave and multi-component migration results is the decoupling of P-and S-waves.The elastic wavefields are separated during the reconstruction of the source and receiver wavefields to obtain pure P-wave and S-wave wavefields [46].The following two methods are employed in this study.

Wave Equation Decoupling Method
Based on the homogeneous isotropic media and the first-order velocity-stress equations, the velocities of separated P-and S-waves in 3D space can be expressed as follows [47]: In Equation ( 3), v px and v sx represent the P-and S-wave components in the x-direction, respectively.Similarly, v py and v sy represent the P-and S-wave components in the ydirection, and v pz and v sz represent the P-and S-wave components in the z-direction, respectively.By solving Equation (3), the wavefields of pure P-and S-waves can be obtained.With this method, the amplitude and phase information of the P-and S-waves can be better retained.

Auxiliary Variables Method
In addition to the six velocity component equations incorporated in the equation for conventional elastic wave, the method of auxiliary variables introduces an additional component to achieve the vector decomposition of P-and S-waves, which is the compressional wave pressure component equation [48,49].The expression of this equation is as follows: In the method of auxiliary variables, the introduction of the τ p component allows for additional wavefield information to be obtained.The τ p component is analogous to the acoustic pressure component in the acoustic wave equation.In the context of elastic waves, it represents the pressure component of the P-wave [50].

Cross-Correlation Imaging Condition
The RTM algorithm, with its significant impact on imaging accuracy, relies heavily on the imaging condition [51,52].Different imaging conditions employed in the imaging process can lead to different imaging results.The following sections introduce four crosscorrelation imaging conditions used in elastic RTM.

Imaging by Vector Velocity Fields
In conventional elastic RTM imaging, the complete wavefield is employed for the imaging process.For multi-component seismic data, the velocity wavefields of the source (v S (x,t)) and the receiver (v R (x,t)) in the subsurface at each position x can be obtained through forward modeling and wavefield extrapolation.The imaging condition utilizing the three-component velocity wavefields can be expressed as follows: where S and R represent the velocity wavefields of the source and receiver, respectively; i and j correspond to the three Cartesian components: x, y, and z; and I ij represents the imaging structure generated by cross-correlating the i-component of the source wavefield with the j-component of the receiver wavefield.T max is the maximum time of the seismic records.This method can introduce crosstalk between P-and S-waves in both the source wavefield and the receiver wavefield.

Imaging by Scalar and Vector Potentials
Under the elastic RTM imaging condition, the second imaging method is implemented.It utilizes the Helmholtz decomposition during the extrapolation of the wavefield, effectively separating the source wavefield and receiver wavefield into pure P-and S-wave wavefields.The expression of the imaging condition is as follows: where S and R represent the velocity wavefields of the source and receiver, respectively; i and j correspond to different wave modes, specifically the P-and S-waves.T max is the maximum time of the seismic records.This method modifies the phase and amplitude characteristics of the original wavefield when employing the Helmholtz decomposition.

Imaging by Pure Vector P-and S-Waves
Using the separation method described in Equation (3), pure P-and S-wave wavefields in various directions can be obtained.Afterward, during the process of wavefield extrapolation, these separated wavefields are utilized as boundary conditions.The imaging results of the P-and S-wave vector components can be acquired through correlating the source wavefield with the receiver wavefield.Figure 1 illustrates the imaging process for this method.
where S and R represent the velocity wavefields of the source and receiver, respectively; i and j correspond to different wave modes, specifically the P-and S-waves.Tmax is the maximum time of the seismic records.This method modifies the phase and amplitude characteristics of the original wavefield when employing the Helmhol decomposition.

Imaging by Pure Vector P-and S-Waves
Using the separation method described in Equation (3), pure P-and S-wave wavefields in various directions can be obtained.Afterward, during the process of wavefield extrapolation, these separated wavefields are utilized as boundary conditions.The imaging results of the P-and S-wave vector components can be acquired through correlating the source wavefield with the receiver wavefield.Figure 1 illustrates the imaging process for this method.The following expression represents the imaging condition for elastic wave RTM using vector P-and S-wave velocity wavefields: where i represents the direction of the source wavefield and receiver wavefield (x, y, z); j represents the wavefield type of the source wavefield and receiver wavefield (P for P-wave, S for S-wave).Tmax is the maximum time of the seismic records.Since the decoupled ve- The following expression represents the imaging condition for elastic wave RTM using vector P-and S-wave velocity wavefields: where i represents the direction of the source wavefield and receiver wavefield (x, y, z); j represents the wavefield type of the source wavefield and receiver wavefield (P for P-wave, S for S-wave).T max is the maximum time of the seismic records.Since the decoupled velocity-stress equations are used in the calculation of P-and S-waves and the input wavefields are either pure P-wave or pure S-wave, the utilization of the imaging condition described in Equation ( 7) leads to minimized interference in the obtained imaging results.

Imaging by Scalarizing Vector Wavefields
Both the pure P-and S-wave wavefields obtained through decoupling are vector fields.However, the results obtained through RTM imaging are not align with the expected reflection coefficients depicted in the imaging profile, lacking clear physical interpretation.To overcome this issue, a scalarization imaging condition is proposed.It converts the vector wavefields obtained from wavefield separation into scalar quantities before crosscorrelation imaging.Figure 2 describes the detailed imaging process of this method.The following expression represents the imaging condition: where PP I , PS I , SP I , and SS I denote PP, PS, SP, and SS images, respectively; P S and S S denote P-and S-wave vector velocity fields at the source side, respectively; and P R and S R denote P-and S-wave vector velocity fields at the receiver side, respectively.

Denoising of the RTM Imaging Results
Although RTM is effective in imaging complex structures, low-frequency artifacts can occur when applying imaging conditions [53,54].On this basis, the Poynting vector and Laplace filtering are proposed to a enuate RTM imaging artifacts.In elastic wavefields, the Poynting vector can be calculated by considering the stress and velocity components of the particles, as follows: The following expression represents the imaging condition: where I PP , I PS , I SP , and I SS denote PP, PS, SP, and SS images, respectively; S P and S S denote P-and S-wave vector velocity fields at the source side, respectively; and R P and R S denote P-and S-wave vector velocity fields at the receiver side, respectively.

Denoising of the RTM Imaging Results
Although RTM is effective in imaging complex structures, low-frequency artifacts can occur when applying imaging conditions [53,54].On this basis, the Poynting vector and Laplace filtering are proposed to attenuate RTM imaging artifacts.In elastic wavefields, the Poynting vector can be calculated by considering the stress and velocity components of the particles, as follows: In the equations, τ ij (i, j = x, y, z) and v i (i = x, y, z) represent the stress and velocity components, respectively.E x , E y , and E z represent the vectors of energy flux density in the x, y, and z directions, respectively.A positive value of E i (i = x, y, z) indicates that the propagation of the wavefield in the positive direction along the corresponding axis.Conversely, a negative value of E i (i = x, y, z) indicates that the propagation of the wavefield in the negative direction along the corresponding axis.In the 3D case, the Laplace operator can be expressed as follows: where ∇ 2 denotes the Laplace operator; k x , k y , and k z represent the wave numbers along the coordinate axes; θ represents the incident angle; V denotes the medium's velocity; and ω is the angular frequency.Following Laplace filtering, the imaging noise is entirely eliminated within the imaging profile when the incident angle is θ = 90 • .When θ < 90 • , the imaging noise is partially suppressed to a certain extent.

Verifications and Discussion
In the following section, we start by discussing the vector decoupling for elastic wave separation.We then validate and discuss the imaging results under four different imaging conditions.Lastly, we verify and discuss the noise suppression methods for RTM imaging.

Discussion on Vector Decoupling for Elastic Wave Separation
A cross-well seismic model was designed to evaluate the effectiveness of the elastic wave separation algorithm in cross-well seismic exploration, as depicted in Figure 3.The effective computational domain of the model is 100.0 m × 100.0 m × 150.0 m.The physical properties of the three-layered medium, arranged from top to bottom, are as follows: the P-wave velocities are 2200.0m/s, 2600.0 m/s, and 3300.0 m/s, respectively; the S-wave velocities are 1200.0m/s, 1450.0 m/s, and 1850.0 m/s, respectively; the densities of the medium are 1800.0kg/m 3 , 2000.0 kg/m 3 , and 2500.0 kg/m 3 , respectively; and the layer thicknesses are 60.0 m, 60.0 m, and 30.0 m, respectively.A sixth-order finite difference algorithm with a spatial grid size of 1.0 m in each dimension was used for the forward modeling simulations.The seismic source employed in the study utilized a Gaussian first derivative wavelet with a dominant frequency of 120.0 Hz.The seismic source type was an explosive source located at coordinates (50.0 m, 50.0 m, 3.0 m).The horizontal position of the receiving well was (20.0 m, 0.0 m), covering the entire well range.The interline spacing between receivers was set at 2.0 m, resulting in a total of 75 receiver points.The recording length was set to 200.0 ms.
The three-component wavefield recordings are depicted in Figure 4. Figure 4a-c represent the v x component, while Figure 4d-f represent the v y component, and Figure 4g-i represent the v z component.Figure 4a illustrates the undecomposed v x component recording, Figure 4d illustrates the undecomposed v y component recording, and Figure 4g illustrates the undecomposed v z component recording.In these wavefields, all three components contain information about both P-and S-waves.However, Figure 4b,e,h display the wavefield recordings of the P-wave in the three components, exclusively capturing P-wave information.Similarly, Figure 4c,f,i show the wavefield recordings of the S-wave in the three components, exclusively containing shear wave information.Due to the use of an explosive seismic source, the recorded direct wave includes a compressional P-wave.Records obtained within the well are more complex than surface records.They encompass transmitted compressional waves TP1 and TP2 through interfaces and reflected compressional waves RPP1 and RPP2 originating from two interfaces.Furthermore, transmitted converted shear waves TPS1 and TPS2 and reflected converted shear waves RPS1 and RPS2 are also obtained, accompanied by interbed multiple conversions.In the wavefield recordings of the P-wave (Figure 4b,e,h), it is evident that the main energy is occupied by P-waves (TP, RPP).However, in the wavefield recordings of the S-wave (Figure 4c,f,i), the energy of P-waves is greatly attenuated, and the S-waves (TPS, RPS) are the predominant component in the wavefield.The wavefield recordings are comprehensive, offering valuable information for understanding the wave propagation process.This observation strongly illustrate that the method employed in this study effectively on vector decoupling for elastic wave separation in the three-component cross-well recorded wavefield.It demonstrates the effectiveness and feasibility of the proposed method.
A single seismic trace was extracted from the inter-well Vz component, Vpz component, and Vsz component at a depth of 30.0 m for comparison.The comparison results are presented in Figure 5.It can be seen that for cross-well seismic recordings, the intermediate variable method successfully decomposes the data into the horizontal P-wave component Vpz, which exclusively contains all the P-wave information of the Vz component.The amplitudes and phase characteristics of the P-wavefields in the Vpz recordings precisely match those of the Vz component.Similarly, the decomposed S-wave component Vsz closely matches the S-wave information in the Vz component.The above results suggest that the method is effective in separating P-and S-waves from the elastic wavefield in a 3D cross-well environment.Notably, amplitude preservation is realized during decomposition, allowing for the isolation of P-and S-waves without interference.This characteristic enables the calculation of cross-well seismic source wavefield propagation directions and improves the accuracy of imaging results during subsequent RTM processes.that the method is effective in separating P-and S-waves from the elastic wave 3D cross-well environment.Notably, amplitude preservation is realized during d sition, allowing for the isolation of P-and S-waves without interference.This ch tic enables the calculation of cross-well seismic source wavefield propagation d and improves the accuracy of imaging results during subsequent RTM processe In order to quantitively analyze the two wavefield separation methods, single-trace records extracted from the Vz component were normalized and compared.The comparison of the single-trace records is illustrated in Figure 6a.It can be seen that the characteristics of both methods are consistent after normalizing the amplitudes.However, the key difference is that the intermediate variable method generates an independent record of the auxiliary variables τ p , as illustrated in Figure 6b.
In order to quantitively analyze the two wavefield separation methods, single-trace records extracted from the Vz component were normalized and compared.The comparison of the single-trace records is illustrated in Figure 6a.It can be seen that the characteristics of both methods are consistent after normalizing the amplitudes.However, the key difference is that the intermediate variable method generates an independent record of the auxiliary variables τp, as illustrated in Figure 6b.

Discussion on the Imaging Results under Various Imaging Conditions
The model depicted in Figure 7

Imaging Results Obtained Utilizing Vector Velocity Fields
Numerical simulation and imaging were conducted using the observational system described above, and the results obtained with the imaging strategy in Section 2.3.1 are shown in Figure 8. Due to the implementation of single-shot recordings, the completeness of the in-well 3D imaging results and the visual depiction of the structure may not be as good as conventional surface-based imaging.However, an interface can be observed at a depth of 60.0 m despite relatively low continuity.This limitation is due to the constraints imposed by the cross-well observation system.I xx and I yy exhibit a significant similarity, differing only by 90 degrees in orientation, which aligns with the expected theoretical results.These two components rely heavily on the v x and v y components, highlighting the imaging results of PS-waves.However, since the seismic source is centered at XOY plane and the receiving well is located at the outer boundary, the central region of the model lacks crucial wavefield information, making it unsuitable for imaging.As a result, imaging can only be performed in the boundary region.The imaging results for I zz are mainly dependent on the v z component, thus primarily representing the imaging results of PP-waves.Similarly, imaging in the central region of the model is not feasible due to the absence of wavefield information.Furthermore, high levels of noise can be observed in the imaging results.This noise originates from the direct correlation imaging of velocity components, leading to interference among different wave modes and a large number of noise artifacts.As a result, the quality of the imaging results is compromised.The lack of a clear physical interpretation of the offset imaging results, represented by the horizontal and vertical components, hinders subsequent tasks such as interpretation and inversion.
sence of wavefield information.Furthermore, high levels of noise can be observed in the imaging results.This noise originates from the direct correlation imaging of velocity components, leading to interference among different wave modes and a large number of noise artifacts.As a result, the quality of the imaging results is compromised.The lack of a clear physical interpretation of the offset imaging results, represented by the horizontal and components, hinders subsequent tasks such as interpretation and inversion.

Imaging Results Obtained Utilizing Scalar and Vector Potentials
Figure 9 illustrates the imaging results obtained using the imaging strategy described in Section 2.3.2.It can be seen that due to the implementation of single-shot recordings and the constraints of the cross-well observation system, the 3D imaging results within the well have relatively poor interface continuity.The imaging results of the PP-waves are close to the source well area, while those of the PS-waves are close to the receiver well area, consistent with the actual situation.Due to the interferences among different mode waves, the imaging results of PS-waves exhibit significantly higher noise levels compared to those of PP-waves.In 3D space, both the imaging interfaces of PP-and PS-waves exhibit polarity reversals.This phenomenon poses challenges for the stacking of multiple shots.

Imaging Results Obtained Utilizing Scalar and Vector Potentials
Figure 9 illustrates the imaging results obtained using the imaging strategy described in Section 2.3.2.It can be seen that due to the implementation of single-shot recordings and the constraints of the cross-well observation system, the 3D imaging results within the well have relatively poor interface continuity.The imaging results of the PP-waves are close to the source well area, while those of the PS-waves are close to the receiver well area, consistent with the actual situation.Due to the interferences among different mode waves, the imaging results of PS-waves exhibit significantly higher noise levels compared to those of PP-waves.In 3D space, both the imaging interfaces of PP-and PS-waves exhibit polarity reversals.This phenomenon poses challenges for the stacking of multiple shots.

Imaging Results Obtained Using Pure Vector P-and S-Wave
The P-wave components (IxP, IyP, and IzP) were imaged in various directions modes using the imaging strategy outlined in Section 2.3.3.The results are presented Figure 10a-c.Additionally, the imaging results for the S-wave components (IxS, IyS, and are depicted in Figure 10d-e.It can be observed that the imaging results of the six com nents in Figure 10 are relatively similar to those in Figure 8.The imaging results of two horizontal components also exhibit a 90° rotation, and the central region of the mo still lacks wavefield information.It can be clearly observed that employing the separa

Imaging Results Obtained Using Pure Vector P-and S-Wave
The P-wave components (I xP , I yP , and I zP ) were imaged in various directions and modes using the imaging strategy outlined in Section 2.3.3.The results are presented in Figure 10a-c.Additionally, the imaging results for the S-wave components (I xS , I yS , and I zS ) are depicted in Figure 10d-e.It can be observed that the imaging results of the six components in Figure 10 are relatively similar to those in Figure 8.The imaging results of the two horizontal components also exhibit a 90 • rotation, and the central region of the model still lacks wavefield information.It can be clearly observed that employing the separated wave equations and utilizing either pure compressional or shear waves as boundary conditions results in images with significantly reduced noise levels compared to Figure 8.The separated wave equations and the input of single-mode waves effectively reduce interferences in imaging, significantly increasing the signal-to-noise ratio of the imaging results.However, different directions and modes produce more images, and the physical significance of the offset imaging results expressed by the two horizontal components and the vertical component is unclear, which is inconducive to the superimposition of the imaging results.

Imaging Results Obtained Using Pure Vector P-and S-Wave
The P-wave components (IxP, IyP, and IzP) were imaged in various directions and modes using the imaging strategy outlined in Section 2.3.3.The results are presented in Figure 10a-c.Additionally, the imaging results for the S-wave components (IxS, IyS, and IzS) are depicted in Figure 10d-e.It can be observed that the imaging results of the six components in Figure 10 are relatively similar to those in Figure 8.The imaging results of the two horizontal components also exhibit a 90° rotation, and the central region of the model still lacks wavefield information.It can be clearly observed that employing the separated wave equations and utilizing either pure compressional or shear waves as boundary conditions results in images with significantly reduced noise levels compared to Figure 8.The separated wave equations and the input of single-mode waves effectively reduce interferences in imaging, significantly increasing the signal-to-noise ratio of the imaging results.However, different directions and modes produce more images, and the physical significance of the offset imaging results expressed by the two horizontal components and the vertical component is unclear, which is inconducive to the superimposition of the imaging results.

Imaging Results Obtained Using Scalarizing Vector Wavefields
In order to achieve easier stacking of multiple shots and clearer physical interpretation of the imaging results, the imaging strategy described in Section 2.3.4 was implemented.This strategy involved converting the vectorial P-waves and vectorial S-waves into scalar quantities.The scalarized imaging results for PP-waves and PS-waves are illustrated in Figure 11.It is evident from Figure 11a,b that the scalarized imaging results of PP-waves and PS-waves are relatively continuous and have no directionality.This characteristic facilitates the stacking of multiple shots after imaging.In addition, the τp component can be extracted from the wavefield using the auxiliary variables method, which is also applicable for imaging purposes.The imaging results obtained using this component are illustrated in Figure 11c.The imaging results for the τp component are similar to the imaging results of PP-waves, with an improved signal-to-noise ratio and continuity.These results can be effectively employed in subsequent applications.

Imaging Results Obtained Using Scalarizing Vector Wavefields
In order to achieve easier stacking of multiple shots and clearer physical interpretation of the imaging results, the imaging strategy described in Section 2.3.4 was implemented.This strategy involved converting the vectorial P-waves and vectorial S-waves into scalar quantities.The scalarized imaging results for PP-waves and PS-waves are illustrated in Figure 11.It is evident from Figure 11a,b that the scalarized imaging results of PP-waves and PS-waves are relatively continuous and have no directionality.This characteristic facilitates the stacking of multiple shots after imaging.In addition, the τ p component can be extracted from the wavefield using the auxiliary variables method, which is also applicable for imaging purposes.The imaging results obtained using this component are illustrated in Figure 11c.The imaging results for the τ p component are similar to the imaging results of PP-waves, with an improved signal-to-noise ratio and continuity.These results can be effectively employed in subsequent applications.
into scalar quantities.The scalarized imaging results for PP-waves and PS-waves are illustrated in Figure 11.It is evident from Figure 11a,b that the scalarized imaging results of PP-waves and PS-waves are relatively continuous and have no directionality.This characteristic facilitates the stacking of multiple shots after imaging.In addition, the τp component can be extracted from the wavefield using the auxiliary variables method, which is also applicable for imaging purposes.The imaging results obtained using this component are illustrated in Figure 11c.The imaging results for the τp component are similar to the imaging results of PP-waves, with an improved signal-to-noise ratio and continuity.These results can be effectively employed in subsequent applications.According to the four different imaging results mentioned above, firstly, the imaging results by vector velocity fields exhibit significant noise, resulting in low-quality images.Furthermore, the imaging results lack clear physical meaning, which hinders subsequent interpretation and inversion tasks.Secondly, the imaging by scalar and vector potentials According to the four different imaging results mentioned above, firstly, the imaging results by vector velocity fields exhibit significant noise, resulting in low-quality images.Furthermore, the imaging results lack clear physical meaning, which hinders subsequent interpretation and inversion tasks.Secondly, the imaging by scalar and vector potentials show relatively poor interface continuity, with polarity reversals observed on both the imaging interfaces of PP waves and PS waves.Moreover, the imaging results by pure vector P-and S-waves exhibit a high-quality representation with minimal noise interference.However, the physical significance of the imaging results is also unclear.Finally, the imaging results by scalarizing vector wavefields exhibit a high degree of continuity and lack directional bias.This characteristic makes them more suitable for superposition, enabling easier overlaying of multiple images and facilitating a clearer physical interpretation of the imaging results.

RTM Imaging Artifacts Attenuation
Figure 12 compares the pre-and post-denoised imaging profiles of single-shot PPwaves.It can be seen that without noise suppression, the imaging profiles exhibit a large amount of low-frequency noise with high energy.This noise is mainly distributed above the interfaces, significantly interfering with the imaging quality of the interfaces.By employing the Poynting vector to suppress the reflected waves and applying Laplace filtering, the low-frequency noise is effectively eliminated.Consequently, the imaging effect is improved, and the actual position and tilt angle of the interface are more accurately reflected.
amount of low-frequency noise with high energy.This noise is mainly distributed abov the interfaces, significantly interfering with the imaging quality of the interfaces.By em ploying the Poynting vector to suppress the reflected waves and applying Laplace filterin the low-frequency noise is effectively eliminated.Consequently, the imaging effect is im proved, and the actual position and tilt angle of the interface are more accurately reflected (a) (b)

Multiple Sensors Cross-Well 3D RTM Image Results and Analysis
In this section, we first discuss the sensor se ings and design of the observation sy tem for 3D cross-well seismic exploration.Subsequently, using the observation system we conduct 3D imaging of layered media model and high-velocity ellipsoidal boulde models.Meanwhile, a comprehensive and detailed analysis of the imaging results is pro vided.

Sensor Se ings and Observation System
The cross-well seismic observation system distinguishes itself from the surface-base 3D seismic methods by enabling the large-scale and high-density deployment of seism acquisition lines.However, the cross-well seismic observation system is influenced by th actual well layout.Based on practical engineering considerations, this paper discusses th layout of the 3D multi-sensors cross-well observation system in commonly encountere square exploration areas, as depicted in Figure 13.Within a 50.0 m × 50.0 m square plan a total of 20 well locations were arranged for the analysis.The specific coordinates of eac well in the XOY plane can be found in Table 1.The data acquisition was performed in th following steps: Initially, the source was subjected to an explosion in Well-1, while th other 19 were used to receive seismic signals.Then, the source was sequentially moved t

Multiple Sensors Cross-Well 3D RTM Image Results and Analysis
In this section, we first discuss the sensor settings and design of the observation system for 3D cross-well seismic exploration.Subsequently, using the observation system, we conduct 3D imaging of layered media model and high-velocity ellipsoidal boulder models.Meanwhile, a comprehensive and detailed analysis of the imaging results is provided.

Sensor Settings and Observation System
The cross-well seismic observation system distinguishes itself from the surface-based 3D seismic methods by enabling the large-scale and high-density deployment of seismic acquisition lines.However, the cross-well seismic observation system is influenced by the actual well layout.Based on practical engineering considerations, this paper discusses the layout of the 3D multi-sensors cross-well observation system in commonly encountered square exploration areas, as depicted in Figure 13.Within a 50.0 m × 50.0 m square plane, a total of 20 well locations were arranged for the analysis.The specific coordinates of each well in the XOY plane can be found in Table 1.The data acquisition was performed in the following steps: Initially, the source was subjected to an explosion in Well-1, while the other 19 were used to receive seismic signals.Then, the source was sequentially moved to the next well, and the remaining 19 wells continued to receive the seismic signals.This process lasted until Well-20 was reached.Table 1.The well locations within the square area in the XOY plane. No.
Well    3 , and an interface depth of 60.0 m.After removing the direct waves from the recordings, the decoupled 3D elastic wave equation RTM can be used to image the layered medium in 3D across walls.The corresponding imaging results are presented in Figure 14b-d.Among them, Figure 14b depicts the imaging results for PP-waves, Figure 14c illustrates the imaging results for PS-waves, and Figure 14d shows the imaging results for τ p .The wavefield energy is precisely concentrated at the actual interface, which is very similar to the real model.This observation result indicates that using cross-well elastic wave data and sensors can ensure accurate imaging of the 3D cross-well layered medium, encompassing both P-and S-waves.The accurate imaging results demonstrate the feasibility of the proposed method.Notably, the source locations are all positioned above the interface, and the imaging results successfully reveal the structural characteristics of the underlying interface.Compared to the traditional first-arrival travel-time imaging method, the RTM imaging method based on reflected waves has a wider exploration range.It eliminates the depth limitation of traditional observing systems for exploration.The physical properties of the ellipsoidal body match those of the second-layer medium.After eliminating the direct waves from the recorded data, the cross-well 3D elastic wave imaging results for the high-velocity ellipsoid boulder can be obtained using the decoupled 3D elastic wave RTM equation.These results are illustrated in Figure 15b-d.Specifically, Figure 15b depicts the imaging result of the PP-wave, Figure 15c illustrates the imaging result of the PS-wave, and Figure 15d shows the imaging result of τ p .It can be seen that the wavefield energy is concentrated at the interfaces and the position of the ellipsoid boulder in the center.The locations of the interfaces align with those of the coherent events, demonstrating a high consistency with the actual model.Based on this result, it can be concluded that the 3D cross-well elastic wave RTM imaging method with multiple wells and sensors can capture the P-and S-wave structures of localized heterogeneous bodies and perform 3D imaging.In addition to horizontal interfaces, the geological structures of lateral media are effectively revealed, and variations in lateral structures are accurately characterized.Moreover, the limitations of traditional 2D imaging methods in imaging lateral media are successfully addressed.In conclusion, the proposed method provides a fundamental basis for extracting dynamic parameters and facilitates the acquisition of comprehensive geological and geophysical information about 3D geological structures.Therefore, this method can be used as a theoretical reference and a solid foundation for detecting the geological formations of high-velocity ellipsoid boulders, especially in the field of engineering exploration.The physical properties of the ellipsoidal body match those of the secondlayer medium.After eliminating the direct waves from the recorded data, the cross-well 3D elastic wave imaging results for the high-velocity ellipsoid boulder can be obtained using the decoupled 3D elastic wave RTM equation.These results are illustrated in Figure 15b-d.Specifically, Figure 15b depicts the imaging result of the PP-wave, Figure 15c illustrates the imaging result of the PS-wave, and Figure 15d shows the imaging result of τp.It can be seen that the wavefield energy is concentrated at the interfaces and the position of ing methods in imaging lateral media are successfully addressed.In conclusion, the proposed method provides a fundamental basis for extracting dynamic parameters and facilitates the acquisition of comprehensive geological and geophysical information about 3D geological structures.Therefore, this method can be used as a theoretical reference and a solid foundation for detecting the geological formations of high-velocity ellipsoid boulders, especially in the field of engineering exploration.

Conclusions
In coastal engineering, a large amount of unsaturated loose sediments (plentiful sand) in the coastal zone leads to a strong a enuation of seismic waves.It poses significant challenges for ground seismic methods to transmit sufficient energy through the weathered zone to reach the desired depth.Compared to surface seismic reflection or refraction surveys, cross-well seismic acquisition methods emerge as an optimal choice for effectively targeting specific coastal areas.However, current studies on cross-well RTM imaging are limited to 2D space or 3D P-wave RTMs.To address the limitations of current studies, this study proposed the 3D cross-well elastic RTM imaging based on a multi-wave and multi-

Conclusions
In coastal engineering, a large amount of unsaturated loose sediments (plentiful sand) in the coastal zone leads to a strong attenuation of seismic waves.It poses significant challenges for ground seismic methods to transmit sufficient energy through the weathered zone to reach the desired depth.Compared to surface seismic reflection or refraction surveys, cross-well seismic acquisition methods emerge as an optimal choice for effectively targeting specific coastal areas.However, current studies on cross-well RTM imaging are limited to 2D space or 3D P-wave RTMs.To address the limitations of current studies, this study proposed the 3D cross-well elastic RTM imaging based on a multi-wave and multi-component technique in coastal engineering exploration.The practical benefit of this method is the utilization of the Earth's elastic medium without increasing costs, which is used to obtain information about the distribution of subsurface rock layers, interfaces, and other structures in coastal engineering projects.Importantly, this can be achieved without resorting to extensive excavation or drilling operations.
Firstly, based on the vector decoupled elastic wave equation, this study compares and analyzes the amplitude-preserving separation algorithms for P-and S-waves implemented using the direct decomposition method and the auxiliary variable method.The wavefield decoupling algorithm achieves accurate separation of vector P-and S-waves without introducing any distortions in amplitude and phase.This feature facilitates accurate imaging of P-and S-waves in subsequent stages.Moreover, the algorithm incorporating intermediate variables integrates the P-wave pressure equation into the traditional elastic wave equations, offering significant advantages for independent imaging in the later stages.Then, this study focuses on analyzing four imaging methods under cross-correlation imaging conditions in RTM.The imaging performances of these four types of imaging strategies are compared and analyzed using experimental models.The results indicate that well-defined P-and S-wave imaging profiles can be obtained using scalarizing vector wavefield imaging conditions, and the imaging results have unique physical significance.Finally, the layer model and high-velocity ellipsoid boulder model were subjected to an RTM imaging experiment.The results indicate that the proposed 3D elastic wave imaging method can effectively generate accurate 3D cross-well profiles of P-and S-waves, thus accurately describing the geological structure.The implementation of multi-wave and multi-component RTM imaging significantly improves the utilization of wavefield information during the imaging process, thereby providing novel insights for cross-well seismic exploration.Moreover, abundant geological and geophysical information can be obtained for analyzing and interpreting 3D geological structures in coastal areas.The findings of this study have crucial theoretical significance and practical implications for exploration and development in real-world production, as well as for the investigation of geological structures in coastal engineering projects.

Figure 1 .
Figure 1.The flowchart illustrating the imaging process using pure vector P-and S-waves.

Figure 1 .
Figure 1.The flowchart illustrating the imaging process using pure vector P-and S-waves.

Figure 2 .
Figure 2. The flowchart illustrating the imaging process using scalarizing vector wavefields.

Figure 2 .
Figure 2. The flowchart illustrating the imaging process using scalarizing vector wavefields.

JFigure 4 .
Figure 4. Cross-well three-component recordings and decomposed P-and S-wave recordings.(a-c) vx component; (d-f) vy component; (g-i) vz component.Figure 4. Cross-well three-component recordings and decomposed P-and S-wave recordings.(a-c) v x component; (d-f) v y component; (g-i) v z component.

Figure 4 .
Figure 4. Cross-well three-component recordings and decomposed P-and S-wave recordings.(a-c) vx component; (d-f) vy component; (g-i) vz component.Figure 4. Cross-well three-component recordings and decomposed P-and S-wave recordings.(a-c) v x component; (d-f) v y component; (g-i) v z component.

Figure 5 . 2 Figure 5 .
Figure 5.The seismic trace comparison between the cross-well Vz component and the V components.

Figure 6 .Figure 6 .
Figure 6.(a) The single-trace comparison plot; (b) the record of the τp component within the well.3.2.Discussion on the Imaging Results under Various Imaging Conditions The model depicted in Figure 7 was established to analyze the validity of each imaging result.It comprises a two-layer model with a computational domain of 100.0 m × 100.0 m × 80.0 m.The physical properties of the two layers from top to bo om are P-wave ve-

Figure 6 .
Figure 6.(a) The single-trace comparison plot; (b) the record of the τp component within th

Figure 7 .
Figure 7. Schematic representation of a layered model and an observational system.The sm rectangle represents the position of the wells.The blue asterisk represents the location of the

Figure 7 .
Figure 7. Schematic representation of a layered model and an observational system.The small red rectangle represents the position of the wells.The blue asterisk represents the location of the source.

Figure 8 .
Figure 8. Imaging results obtained using vector velocity fields for all three components.(a) Ixx component imaging results; (b) Iyy component imaging results; (c) Izz component imaging results.

Figure 8 .
Figure 8. Imaging results obtained using vector velocity fields for all three components.(a) I xx component imaging results; (b) I yy component imaging results; (c) I zz component imaging results.

Figure 9 .
Figure 9. Imaging results obtained utilizing scalar and vector potentials.(a) PP-wave imaging results; (b) PS-wave imaging results.

Figure 10 .
Figure 10.Imaging results obtained using pure vector P-and S-waves.(a) I xP component imaging results; (b) I yP component imaging results; (c) I zP component imaging results; (d) I xS component imaging results; (e) I yS component imaging results; (f) I zS component imaging results.

Figure 11 .
Figure 11.Imaging results obtained utilizing pure vector P-and S-waves.(a) PP-wave imaging results; (b) PS-wave imaging results; (c) τp component imaging results.

Figure 11 .
Figure 11.Imaging results obtained utilizing pure vector P-and S-waves.(a) PP-wave imaging results; (b) PS-wave imaging results; (c) τp component imaging results.

Figure 12 .
Figure 12.Pre-and post-denoised imaging profiles of single-shot PP-waves: (a) image results before denoising; (b) image results after denoising.
J. Mar.Sci.Eng.2024, 12, x FOR PEER REVIEW 16 of 22 the next well, and the remaining 19 wells continued to receive the seismic signals.This process lasted until Well-20 was reached.

Figure 13 .
Figure 13.Three-dimensional multi-cross-well observation system diagram.The small red rectangle represents the position of the wells.

Figure 13 .
Figure 13.Three-dimensional multi-cross-well observation system diagram.The small red rectangle represents the position of the wells.

2 .
Figure14aillustrates a 3D layered medium model.The physical parameters of the two layers from top to bottom are as follows: P-wave velocities of 2200.0 m/s and 2600.0 m/s, Swave velocities of 1200.0 m/s and 1450.0 m/s, densities of 1800.0 kg/m 3 and 2000.0kg/m 3 , and an interface depth of 60.0 m.After removing the direct waves from the recordings, the decoupled 3D elastic wave equation RTM can be used to image the layered medium in 3D across walls.The corresponding imaging results are presented in Figure14b-d.Among them, Figure14bdepicts the imaging results for PP-waves, Figure14cillustrates the imaging results for PS-waves, and Figure14dshows the imaging results for τ p .The wavefield energy is precisely concentrated at the actual interface, which is very similar to the real model.This observation result indicates that using cross-well elastic wave data and sensors can ensure accurate imaging of the 3D cross-well layered medium, encompassing both P-and S-waves.The accurate imaging results demonstrate the feasibility of the proposed method.Notably, the source locations are all positioned above the interface, and the imaging results successfully reveal the structural characteristics of the underlying interface.Compared to the traditional first-arrival travel-time imaging method, the RTM imaging method based on reflected waves has a wider exploration range.It eliminates the depth limitation of traditional observing systems for exploration.

Figure 15
Figure 15 displays the high-velocity ellipsoid boulder model.The physical parameters of the two layers from top to bottom are as follows: P-wave velocities of 2200.0 m/s and 2600.0 m/s, S-wave velocities of 1200.0 m/s and 1450.0 m/s, densities of 1800.0 kg/m 3 and 2000.0kg/m 3 , and an interface depth of 60.0 m.The high-speed ellipsoidal boulder has three axis radii: a = b = 5.0 m, c = 10.0 m.The ellipsoidal body is centered at (25.0 m, 25.0 m, 30.0 m).The physical properties of the ellipsoidal body match those of the second-layer medium.After eliminating the direct waves from the recorded data, the cross-well 3D elastic wave imaging results for the high-velocity ellipsoid boulder can be obtained using the decoupled 3D elastic wave RTM equation.These results are illustrated in Figure15b-d.Specifically, Figure15bdepicts the imaging result of the PP-wave, Figure15cillustrates the imaging result of the PS-wave, and Figure15dshows the imaging result of τ p .It can be seen that the wavefield energy is concentrated at the interfaces and the position of the ellipsoid boulder in the center.The locations of the interfaces align with those of the coherent events, demonstrating a high consistency with the actual model.Based on this result, it can be concluded that the 3D cross-well elastic wave RTM imaging method with multiple wells and sensors can capture the P-and S-wave structures of localized heterogeneous bodies and perform 3D imaging.In addition to horizontal interfaces, the geological structures of lateral media are effectively revealed, and variations in lateral structures are accurately characterized.Moreover, the limitations of traditional 2D imaging methods in imaging lateral media are successfully addressed.In conclusion, the proposed method provides a fundamental basis for extracting dynamic parameters and facilitates the acquisition of comprehensive geological and geophysical information about 3D geological structures.Therefore, this method can be used as a theoretical reference and a solid foundation for detecting the geological formations of high-velocity ellipsoid boulders, especially in the field of engineering exploration.

Figure 14 .
Figure 14.Three-dimensional cross-well layered medium model and the elastic waves RTM imaging results.(a) Layered medium model; (b) PP-waves imaging results; (c) PS-waves imaging results; (d) τp component imaging results.

4. 3 .
Figure 15 displays the high-velocity ellipsoid boulder model.The physical parameters of the two layers from top to bo om are as follows: P-wave velocities of 2200.0 m/s and 2600.0 m/s, S-wave velocities of 1200.0 m/s and 1450.0 m/s, densities of 1800.0 kg/m 3 and 2000.0kg/m 3 , and an interface depth of 60.0 m.The high-speed ellipsoidal boulder has three axis radii: a = b = 5.0 m, c = 10.0 m.The ellipsoidal body is centered at (25.0 m, 25.0 m, 30.0 m).The physical properties of the ellipsoidal body match those of the secondlayer medium.After eliminating the direct waves from the recorded data, the cross-well 3D elastic wave imaging results for the high-velocity ellipsoid boulder can be obtained using the decoupled 3D elastic wave RTM equation.These results are illustrated in Figure15b-d.Specifically, Figure15bdepicts the imaging result of the PP-wave, Figure15cillustrates the imaging result of the PS-wave, and Figure15dshows the imaging result of τp.It can be seen that the wavefield energy is concentrated at the interfaces and the position of

Figure 14 .
Figure 14.Three-dimensional cross-well layered medium model and the elastic waves RTM imaging results.(a) Layered medium model; (b) PP-waves imaging results; (c) PS-waves imaging results; (d) τp component imaging results.

Figure 15 .
Figure 15.Three-dimensional cross-well high-velocity ellipsoid boulder model and the elastic waves RTM imaging results.(a) High-velocity ellipsoid boulder model; (b) PP-waves imaging results; (c) PS-waves imaging results; (d) τp component imaging results.

Figure 15 .
Figure 15.Three-dimensional cross-well high-velocity ellipsoid boulder model and the elastic waves RTM imaging results.(a) High-velocity ellipsoid boulder model; (b) PP-waves imaging results; (c) PS-waves imaging results; (d) τp component imaging results.

Table 1 .
The well locations within the square area in the XOY plane.