Forward Wavefield Research Articles

Full-Waveform Inversion Imaging of Cortical Bone Using Phased Array Tomography.

Classic ultrasound bone imaging modalities usually demand either a prior knowledge or an advanced estimation on speed of sound (SoS), which not only renders to a burdensome imaging process but also supplies a limited resolution. To overcome these drawbacks, this article proposed a frequency-domain full-waveform inversion (FDFWI) modality using phased array tomography for high-accuracy cortical bone imaging. A transmission scenario of ultrasound wave in 2-D space was presented in the frequency domain to simulate the forward wavefield propagation. Iterations in the inversion process were performed by matching the simulation wavefield to the experimental one from low to high discrete frequency points. Moreover, the association between the maximum initial frequency and the initial SoS model was explored to prevent the occurrence of cycle-skipping phenomenon, which could lead to the outcomes being trapped in local minima. The feasibility and effectiveness of the proposed imaging scheme were testified by simulation, phantom, and ex-vivo studies, with mean relative errors of cortical part being 3.18%, 8.71%, and 9.36%, respectively. It is verified that the proposed FDFWI method is an effective way for parametric imaging of cortical bone without any prior knowledge of sound speed.

Full-waveform ultrasonic imaging for bone

Abstract This paper suggests a full-waveform ultrasonic imaging strategy for bone. A novel forward wavefield propagation is performed in the frequency-domain. To improve the efficiency of full-waveform inversion bone imaging, all relative computations are also conducted in the frequency-domain. Additionally, a conjugate gradient algorithm is employed in the inversion procedure to accelerate the reconstruction of the bone acoustic velocity. Eventually, numerical experiments on different bone models are carried out to demonstrate the effectiveness of explored imaging strategy.

Waveform inversion with structural regularizing constraint based on gradient decomposition

Abstract Full waveform inversion (FWI) can simultaneously update low-to-medium wavenumber velocity components and high-wavenumber velocity components. However, if seismic data lack large-offset data and effective low-frequency components, FWI updates will be dominated by high-wavenumber velocity perturbation. Meanwhile, providing that the initial model is inaccurate, inversion will have the problem of local minima. In this study, FWI is developed with structural regularizing constraint based on gradient decomposition (RGDFWI). By correlating the separated forward wavefield and backward wavefield with specific propagating direction, FWI gradient is decomposed into tomography-mode gradient and migration-mode gradient. We propose an optimized strategy taking full advantage of the two modes of FWI gradient. On the one hand, we use tomography-mode gradient to enhance low-to-medium wavenumber updates. On the other hand, we use migration-mode gradient to apply structural regularizing constraint by estimating structure dip and adding sparsity constraint in Seislet domain. During the inversion process, high-wavenumber structural information constrains and guides low-wavenumber model updates. The results of two numerical tests, Marmousi model test and Overthrust model test, validate the optimized strategy, which can produce a better initial velocity model for FWI. The inversion finally generates a high-precision and high-resolution velocity model.

Improving full-waveform inversion based on sparse regularization for geophysical data

Abstract Full-waveform inversion (FWI) is an advanced geophysical inversion technique. FWI provides images of subsurface structures with higher resolution in fields such as oil exploration and geology. The conventional algorithm minimizes the misfit error by calculating the least squares of the wavefield solutions between observed data and simulated data, followed by gradient direction and model update increment. Since the gradient is calculated by forward and backward wavefields, the high-accuracy model update relies on accurate forward and backward wavefield modelling. However, the quality of wavefield solutions obtained in practical situations could be poor and does not meet the requirements of high-resolution FWI. Specifically, the low-frequency wavefield is easily affected by noise and downsampling, which influences data quality, whereas the high-frequency wavefield is susceptible to spatial aliasing effects that produce imaging artefacts. Therefore, we propose using an algorithm called sparse relaxation regularized regression to optimize the wavefield solution in frequency-domain FWI, which is the forward and backward wavefield obtained from the Helmholtz equation, thus improving FWI's accuracy. The sparse relaxation regularized regression algorithm combines sparsity and regularization, allowing the broadband FWI to reduce the effects of noise and outliers, which can provide data supplementation in the low-frequency band and anti-aliasing in the high-frequency band. Our numerical examples demonstrate the wavefield optimization effect of the sparse relaxation regularized regression-based algorithm in various cases. The improved algorithm's accuracy and stability are verified compared to the Tikhonov regularization algorithm.

Excitation time imaging condition reverse time migration based on physics-informed neural network traveltime calculation with wavefield decomposition using optical flow vector

Abstract Although the excitation-time imaging condition offers a lower memory consumption and higher computational efficiency compared to cross-correlation imaging condition, it has not been widely used in industrial applications because of the accuracy problem of traveltime calculation and the influence of low-wave-number noise. In this paper, we introduce the physics-informed neural network (PINN) algorithm to achieve a high-precision traveltime calculation of the source forward wavefield. Subsequently, we introduce a technique for high-precision wavefield decomposition of the reverse-time wavefield via the optical flow vector, enabling us to realize a correlation-weighted stacking imaging of each wavefield. Model experiments and real data processing show that the proposed traveltime calculation algorithm based on PINN offers high accuracy and good applicability in the excitation time reverse-time migration imaging of complex models, and correlation-weighted stacking imaging based on optical flow vector-based wavefield separation can significantly suppress the noise with low wavenumber and achieve high-precision imaging of complex models.

Crosswell frequency-domain reverse time migration imaging with wavefield decomposition

Abstract Crosswell seismic technique can provide high-resolution imaging and monitoring of the subsurface. However, compared to the surface seismic, crosswell data contain more complex wave components, which increases the difficulty of seismic processing and migrations. Conventional acoustic reverse time migration (RTM) mainly uses the cross-correlation of the source forward and the receiver backward wavefields. The redundant information generated by cross-correlation may undermine the imaging reliability in the crosswell models. Thus, we develop a novel wavefield decomposition imaging condition and only used cross-correlation information of incident and reflection waves in the same propagation directions, which eliminates the artifacts generated from the cross-correlation of wave information unrelated to reflections in the crosswell image. We perform RTM in the frequency domain to maintain efficiency in multi-shot problems. A mono-frequency wavefield decomposition method is applied to separate and process the seismic data. The forward and backward wavefields are reclassified into the up- and downpropagating components. The L1 norm is introduced to enhance the robustness of the proposed imaging method. We then use this method to synthesize data from layering models and analyze the imaging results generated from each pair of cross-correlations using source and receiver wavefields. Results show that the cross-correlation information belonging to the same propagation contributes most to the crosswell image, and the other information always generates migration noises. Moreover, we apply the proposed method to a real-field dataset. Processing results validate the effectiveness of the proposed means for eliminating false events in the crosswell models and improving image quality.

Viscoacoustic reverse time migration with stable and effective two-way attenuation compensation

Intrinsic attenuation in seismic wave propagation leads to amplitude dissipation and phase dispersion in seismic data. The attenuation (∼1/ Q) effect is a common cause for inaccurate image locations and dimmed image energy in reverse time migration (RTM). Conventional viscoacoustic or viscoelastic RTM methods implement attenuation compensation by reversing the sign of the amplitude loss term and keeping the dispersion term unchanged for forward and backward wavefield extrapolation, which requires the decoupled viscoacoustic wave equation and gives rise to the numerical instability issue. To address these problems, we develop a stable and effective two-way attenuation-compensated viscoacoustic RTM method. Our method uses a new two-way normalized crosscorrelation imaging condition to compensate the attenuation effect. In the new imaging condition, only the attenuated forward wavefield of sources and the viscoacoustic Green’s functions of receivers are required, which eliminates the instability source arising from the operation of reversing the sign of the dissipation term in the viscoacoustic wave equation. Because no modifications are made to the viscoacoustic wave equation, the new imaging condition can be calculated without any additional computational cost and can be applied to decoupled and coupled viscoacoustic wave equation-based RTMs, which is more flexible than other attenuation compensation strategies. Two synthetic tests and one field data example are presented to demonstrate the stability and effectiveness of our method.

Isogeometric multi-resolution full waveform inversion based on the finite cell method

Full waveform inversion (FWI) is an iterative identification process that serves to minimize the misfit of model-based simulated and experimentally measured wave field data. Its goal is to identify a field of parameters for a given physical object. For many years, FWI is very successful in seismic imaging to deduce velocity models of the earth or of local geophysical exploration areas. FWI has also been successfully applied in various other fields, including non-destructive testing (NDT) and biomedical imaging. The inverse optimization process of FWI relies on forward and backward solutions of the (elastic or acoustic) wave equation, as well as on efficient computations of adequate optimization directions. Many approaches employ (low order) finite element or finite difference methods, using parameterized material fields whose resolution is chosen in relation to the elements or nodes of the discretized wave field. In our previous paper (Bürchner et al., 2023), we investigated the potential of using the finite cell method (FCM) as the wave field solver. The FCM offers the advantage that highly complex geometric models can be incorporated easily. Furthermore, we demonstrated that the identification of the model’s density outperforms that of the velocity — particularly in cases where unknown voids characterized by homogeneous Neumann boundary conditions need to be detected. The paper at hand extends this previous study in the following aspects: The isogeometric finite cell analysis (IGA-FCM) – a combination of isogeometric analysis (IGA) and FCM – is applied as the wave field solver, with the advantage that the polynomial degree and subsequently also the sampling frequency of the wave field can be increased quite easily. Since the inversion efficiency strongly depends on the accuracy of the forward and backward wave field solutions and of the gradients of the functional, consistent, and lumped mass matrix discretization are compared. The resolution of the grid describing the unknown material density – thus allowing to identify voids in a physical object – is then decoupled from the knot span grid. Finally, we propose an adaptive multi-resolution algorithm that locally refines the material grid using an image processing-based refinement indicator. The developed inversion framework allows fast and memory-efficient wave simulations and object identification. While we study the general behavior of the proposed approach using 2D benchmark problems, a final 3D problem shows that it can also be used to identify void regions in geometrically complex spatial structures.

Q-Compensated Gaussian Beam Migration under the Condition of Irregular Surface

The viscosity of actual underground media can cause amplitude attenuation and phase distortion of seismic waves. When seismic images are processed assuming elastic media, the imaging accuracy for the deep reflective layer is often reduced. If this attenuation effect is compensated, the imaging quality of the seismic data can be significantly improved. Q-compensated Gaussian beam migration (Q-GBM) is an effective seismic imaging method for viscous media, and it has the advantages of both wave equation and ray-based Q-compensated imaging methods. This study develops a Q-GBM method in visco-acoustic media with an irregular surface. Initially, the basic principles of Gaussian beam in visco-acoustic media are introduced. Then, by correcting the complex-value time of the Gaussian beam in visco-acoustic media, energy compensation and phase correction are carried out for the forward continuation wavefield at the seismic source of the irregular surface and the reverse continuation wavefield at the beam center, which effectively compensates the absorption and attenuation effects of visco-acoustic media on the seismic wavefield. Further, a Q-GBM method under the irregular surface is proposed using cross-correlation imaging conditions. Through migration tests for three numerical models of visco-acoustic media with irregular surfaces, it is verified that our method is an effective depth domain imaging technique for seismic data in visco-acoustic media under the condition of irregular surfaces.

Compensating attenuation effects in full-waveform inversion with dissipation-dispersion decoupling

Seismic attenuation poses challenges to velocity model building from acoustic/elastic full-waveform inversion (FWI). In particular, when constructing the FWI gradient for velocity inversion, recent studies have indicated that the high-attenuation structure could distort the gradient by damping the amplitudes and shifting the kinematic phases, resulting in an imbalanced update and thus unreliable velocity model. These Q effects are particularly significant in the reflection acquisition geometry due to the “double-damping” issue. Here, we develop a Q-compensated FWI algorithm for constructing a Q-free FWI gradient. By using a recently developed viscoacoustic wave propagator, this compensation can be done conveniently by keeping the dispersion term and flipping the dissipation term in the wave equation when we simulate the forward and time-reversed adjoint wavefields. The resultant gradient obtained by interacting these two wavefields has correct kinematics and Q-free amplitudes. This Q-compensated FWI can balance the update between Q- and no- Q-areas, which we determine using synthetic examples. In addition, we illustrate how to take advantage of the dissipation-dispersion decoupling to determine the anomalous Q value in the Q model building workflow via Q-compensated reverse time migration.

Least‐squares reverse time migration using the most energetic source wavefields based on excitation amplitude imaging condition

AbstractLeast‐squares reverse time migration, a linearized inversion problem, can provide high‐quality migration image by minimizing the misfit function, which is defined by predicted and observed data. According to the theory of data‐domain least‐squares reverse time migration, a forward source wavefield that is simulated with a fixed background velocity does not change during iterations. However, storing the forward source wavefield directly into computer memory involves substantial memory consumption. Although a source wavefield reconstruction technique can be applied during least‐squares reverse time migration iterations, this approach can increase the computational cost because of the need for additional wavefield simulations. To alleviate this computational issue in storing the forward source wavefield, we propose an efficient least‐squares reverse time migration scheme based on an excitation amplitude method. Unlike conventional excitation amplitude imaging conditions, the proposed least‐squares reverse time migration scheme enables one to reconstruct the forward source wavefield by convolving a source wavelet with the excitation amplitude of Green's function at the excitation time. With this excitation amplitude method, the forward source wavefield can be efficiently stored in the computer memory because the sizes of the excitation amplitude and excitation time maps are equal to the size of one snapshot. To validate the feasibility of our least‐squares reverse time migration scheme, we perform dot‐product tests and compare forward source wavefields, demigrated data and gradient vectors obtained by conventional least‐squares reverse time migration and our proposed least‐squares reverse time migration. Using numerical tests on synthetic data, we confirm that our least‐squares reverse time migration can produce high‐quality migration results with a significant improvement in computational efficiency with respect to performance time and memory consumption.

Deep crustal structures with reverse time migration applied to offshore wide-angle seismic data: Equatorial and North-West Brazilian margins

We image the Moho discontinuity and deep crustal layers by applying the Reverse Time Migration (RTM) method to wide-angle seismic (WAS) data that were acquired along two different profiles in NW and SE offshore Brazil, where the spacing between OBS is 12.5 km and 150 m between shots. The application of this method is quite uncommon to ocean bottom seismometers (OBS) data due to the OBS wide spacing deployment and low folds. We analyze the effectiveness of the RTM method when applied to the reflectivity of the WAS data. We imaged the structures by cross-correlating the forward and backward wavefields given by the acoustic seismic equation. This allowed us to use each interface crossed by a ray as both a source and a receiver. The velocity models used to perform RTM were previously obtained by applying a procedure of two-dimensional forward ray-tracing followed by a damped least-squares travel time inversion. The results obtained have an unexpected large contribution from the wavefield traveling as refractions within the Earth and, because of that, we can talk about recovering refractors from the RTM results. We obtained strong and continuous refractors for depths of 7–15 km that correspond to the basement and the Moho discontinuity. As we move landwards, the refractors that correspond to the Moho discontinuity disappear due to the deepening of the surface. Finally, we investigate how the slight change on the velocity model influences the result of the RTM method. The obtained results are promising for a wide range of applications at a crustal scale seismic exploration, with wide-angle seismic data.

Modeling and inversion in acoustic-elastic coupled media using energy-stable summation-by-parts operators

Seismic data acquired at the seafloor are valuable in characterizing the subsurface and monitoring producing hydrocarbon fields. To fully use such data, a stable, accurate, and efficient numerical scheme is needed that accounts for acoustic and elastic wave propagation, their interaction at the seafloor interface, and for sources and receivers placed on either side of that interface. Existing methods either make incorrect assumptions or have high implementation and computational costs. We have developed a high-order finite-difference summation-by-parts framework for the acoustic-elastic wave equations in second-order form (in terms of displacements in the solid and an extension of velocity potential in the fluid). Our modified discretization of the elastic operator overcomes the dispersion errors known to plague displacement-based schemes in the high VP/ VS limit. We weakly impose boundary and interface conditions using simultaneous approximation terms, leading to an energy-stable numerical scheme that rigorously handles point injection and extraction. The fully discrete system is self-adjoint after a time reversal and a sign flip and is furthermore a high-order accurate discretization of the continuous problem. The self-adjointness ensures that forward and adjoint wavefields are computed with similar accuracy and simplifies gradient computation for inversion purposes. We find that our numerical scheme achieves accuracy comparable to a more computationally expensive spectral-element method and demonstrate its application to full-waveform inversion using the Marmousi2 model.

3C-3D tunnel seismic reverse time migration imaging: A case study of Pearl River Delta Water Resources Allocation Project

In recent years, tunnel boring machines (TBMs) have been widely used because of their fast construction speed and economical benefits, and the need to obtain geological information ahead of the tunnel to ensure the TBM tunnel construction safety has also arisen. Faced with the complex geological environment, seismic forward-prospecting is the primary method for predicting the unfavorable geological structure in front of the tunnel face. As a crucial step in seismic forward-prospecting method, migration can image geological characteristics and help learn the distribution of unfavorable geology. To achieve all-round imaging of unfavorable geology, we suggested a three-component (3C) three-dimensional (3D) tunnel seismic reverse time migration imaging method. The acoustic wave equation is employed for the forward and backward wavefield extrapolation after the P/S-wave separation to eliminate the imaging artifacts from severe converted waves in tunnels. GPU parallelization is used to accelerate the calculation to further adapt to the rapid excavation of TBMs. The proposed method is validated in two typical numerical models, especially for 3C seismic data that provides all-round geological information to achieve enhanced imaging energy and artifact reduction, which demonstrates the significance of 3C seismic data to migration imaging. This case study was conducted in the Pearl River Delta Water Resources Allocation Project in Guangdong Province, China. The reliability of 3C3D tunnel seismic reverse time migration imaging in field tunnel engineering was investigated. Two consecutive detections were performed, and then the spatial positions of faults and fractures in front of the tunnel face were accurately judged by the combination of the two imaging results. Based on the detection results, the excavation rate was slowed down and the tunnel lining support was increased. The TBM then successfully passed through the risk areas. This case study can serve as a reference for other tunnels with similar concerns regarding the necessity of 3C seismic data and consecutive detections.

Generating surface-offset common-image gathers with backward wavefield synthesis

Surface-offset common-image gathers (CIGs) are an important data format for seismic velocity analysis. However, the reverse time migration (RTM) method, which is wavefield propagation based, does not directly produce surface-offset CIGs because it propagates waves from all the offsets together. Here, we implement the generation of surface-offset CIGs by synthesizing the backward wavefields using the forward wavefields at locations where source and receiver locations overlap. This method produces surface-offset CIGs with high accuracy and low cost when compared with those generated by backpropagating each trace separately. We adopt nonnegative least-squares filters for sparse linear deconvolution. The computational effectiveness of the proposed method increases for higher dimensions, higher-order stencils, and more complicated wave equations. The proposed method works stably on a realistic towed-streamer acquisition system with moderate geometric positioning errors between the locations of airguns and hydrophones.

Attenuation compensation for wavefield‐separation‐based least‐squares reverse time migration in viscoelastic media

ABSTRACTElastic least‐squares reverse time migration can image the multicomponent seismic recordings. However, on the one hand, without considering the intrinsic attenuation of the subsurface, it may produce blurred reflectivity images with incorrect positions of reflectors for seismic recordings with strong attenuation. On the other hand, the crosstalk artefacts created by the different wave modes may severely degrade the imaging quality. To alleviate these crosstalk artefacts and compensate for the attenuation effects, we have developed a time‐domain wavefield‐separation‐based least‐squares reverse time migration approach in viscoelastic media, which utilizes the viscoelastic wave equation on the basis of the standard linear solid model to simulate intrinsic subsurface attenuation. The key point of the proposed approach is that we utilize the separated gradient contribution of PP‐ and PS‐wave modes based on the wavefield separation technique in viscoelastic media to construct the P‐wave and S‐wave velocity images, respectively. Unlike the conventional wavefield‐separation‐based least‐squares reverse time migration approach, which generally uses the new stress–velocity equations to formulate, our proposed wavefield separation scheme fully depends on the conventional viscoelastic wave equation and its adjoint wave equation. By introducing the pure P‐wave stress in the forward and adjoint wavefields, the coupled wavefields can be well decomposed, which requires much less computational costs than the wavenumber‐domain wavefield separation scheme. Numerical examples using the layered and Marmousi‐2 models have shown that the proposed approach can improve the migration image quality under geological conditions with strong attenuation where elastic least‐squares reverse time migration may produce blurred and unfocused events. Meanwhile, the proposed least‐squares reverse time migration approach with the wavefield separation scheme has a better convergence rate and produces fewer crosstalk artefacts than that without the wavefield separation scheme.

Feature Extraction of Music Signal Based on Adaptive Wave Equation Inversion

The digitization, analysis, and processing technology of music signals are the core of digital music technology. There is generally a preprocessing process before the music signal processing. The preprocessing process usually includes antialiasing filtering, digitization, preemphasis, windowing, and framing. Songs in the popular wav format and MP3 format on the Internet are all songs that have been processed by digital technology and do not need to be digitalized. Preprocessing can affect the effectiveness and reliability of the feature parameter extraction of music signals. Since the music signal is a kind of voice signal, the processing of the voice is also applicable to the music signal. In the study of adaptive wave equation inversion, the traditional full-wave equation inversion uses the minimum mean square error between real data and simulated data as the objective function. The gradient direction is determined by the cross-correlation of the back propagation residual wave field and the forward simulation wave field with respect to the second derivative of time. When there is a big gap between the initial model and the formal model, the phenomenon of cycle jumping will inevitably appear. In this paper, adaptive wave equation inversion is used. This method adopts the idea of penalty function and introduces the Wiener filter to establish a dual objective function for the phase difference that appears in the inversion. This article discusses the calculation formulas of the accompanying source, gradient, and iteration step length and uses the conjugate gradient method to iteratively reduce the phase difference. In the test function group and the recorded music signal library, a large number of simulation experiments and comparative analysis of the music signal recognition experiment were performed on the extracted features, which verified the time-frequency analysis performance of the wave equation inversion and the improvement of the decomposition algorithm. The features extracted by the wave equation inversion have a higher recognition rate than the features extracted based on the standard decomposition algorithm, which verifies that the wave equation inversion has a better decomposition ability.

Correlative least-squares reverse time migration in viscoelastic media

Standard elastic least-squares reverse time migration (LSRTM) generally ignores the subsurface attenuation effects which can lead to blurred and unfocused migration images. To correct these attenuation effects, we develop a correlative Q-LSRTM approach in viscoelastic media. We use the viscoelastic wave equation based on the standard linear solid model to simulate the subsurface attenuation. The Born modeling and adjoint state equations in viscoelastic media can be derived by the Born approximation and adjoint state method, respectively. Instead of the conventional amplitude-matching misfit function, the proposed Q-LSRTM approach utilizes the normalized zero-lag crosscorrelation misfit function to emphasize the phase measure between the predicted and observed data. The gradients of the crosscorrelation misfit function with regard to reflectivity images of P- and S-wave impedance are estimated by a zero-lag crosscorrelation between the forward background wavefields and the reversal adjoint wavefields. In view of computational efficiency, we also derive an analytical step-length formula for the correlative Q-LSRTM approach, which only requires to solve the Born modeling operator once. Numerical examples using several 2D models have determined that the proposed Q-LSRTM can improve better image quality, compared with the standard elastic LSRTM approach. Furthermore, the crosscorrelation misfit function is less sensitive to amplitude errors.

A new decoupling and elastic propagator for efficient elastic reverse time migration

Methods to decompose the elastic wavefield into compressional wave (P-wave) and shear wave (S-wave) components in heterogeneous media without wavefield distortions or energy leakage are the key issues in elastic imaging and inversion. We have introduced a decoupled P- and S-wave propagator to form an efficient elastic reverse time migration (RTM) framework, without assuming homogeneous Lamé parameters. Also, no wave-mode conversions occur using the proposed propagator in the presence of strong heterogeneities, which avoids the potential imaging artifacts caused by wave-mode conversions in the receiver-side backward extrapolation. In the proposed elastic RTM framework, the source-side forward wavefield is simulated with a P-wave propagator. The receiver-side wavefield is back extrapolated with the proposed propagator, using the recorded multicomponent seismic data as input. Compared to the conventional elastic RTM, the proposed framework reduces the computational complexity while preserving the imaging accuracy. We have determined its accuracy and efficiency using two synthetic examples.

A reverse time migration-based multistep angular spectrum approach for ultrasonic imaging of specimens with irregular surfaces

We develop a new ultrasonic imaging framework for non-destructive testing of an immersed specimen featuring an irregular top surface and demonstrate its capability of accurately depicting the lower surfaces of multiple damages hidden in the specimen. Central to the framework is a multistep angular spectrum approach (ASA), via which the forward propagation wavefields of wave sources and backward propagation wavefields of the received wave signals are calculated. Upon applying a zero-lag cross-correlation imaging condition of reverse time migration (RTM) to the obtained forward and backward wavefields, the image of the specimen with an irregular surface can be reconstructed, in which hidden damages, if any and regardless of quantity, are visualized. The effectiveness and accuracy of the framework are examined using numerical simulation, followed with experiment, in both of which multiple side-drilled holes, at different locations in aluminum blocks with various irregular surfaces, are characterized. Results have proven that multiple damages in a specimen with an irregular surface can be individually localized, and the lower surface of each damage can further be imaged accurately, thanks to the RTM-based algorithm in which multiple wave reflections from the specimen bottom are taken into wavefield extrapolation. The proposed imaging approach presents higher computational efficiency, compared to conventional RTM, and enhanced imaging contrast over prevailing total focusing methods.