Multi-scale ultrasonic imaging of sub-surface concrete defects

This paper introduces a new method for creating the ultrasonic image of concrete that utilizes both linear and nonlinear wave properties in a single measurement. Linear ultrasonic imaging relies on wave velocity. The resolution of linear ultrasonics for various inclusion sizes and inclusion-to-base material ratios is numerically and experimentally studied. The nonlinear ultrasonic testing (NLUT) image is based on the acoustic nonlinearity coefficient, calculated using inputs of the first and second harmonic amplitudes, wave number, and the distance between the transmitter and receiver. The ultrasonic waveform is decomposed into its harmonics using wavelet transform. Their wave velocities, wave numbers, and amplitudes are extracted to calculate the linear wave velocity and the acoustic nonlinearity coefficient. Four calibration samples with different inclusions are tested to evaluate the effectiveness of combining two imaging methods: pure concrete, concrete with small foam inserts simulating air inclusions, and two concrete samples with larger foam inclusions. The angle dependence of acoustic nonlinearity coefficient is shown. The results show that constructing both linear and nonlinear ultrasonic imaging is necessary when the defect’s properties are unknown.


Introduction
Concrete is one of the oldest construction materials used to build various structural systems such as roads, bridges, tunnels, dams, and power plants [1].Concrete structures may develop surface and subsurface defects due to environmental factors, extreme events, overload, or aging.Surface defects can be detected with visual inspection or image processing methods using thermal or infrared images.On the other hand, detecting subsurface defects such as corrosion, porosity, and Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.internal concrete crashing due to excessive shear is still a challenging problem.Linear and nonlinear ultrasonics have been widely adapted to detect defects in concrete.Defect detection in concrete using linear ultrasonics is based on the change in attenuation [2], wave velocity or time of flight (TOF) [3][4][5], and wave modality [6] concerning concrete conditions.
While most ultrasonic methods applied to concrete involve local measurements, generating ultrasonic images of concrete is crucial for localizing defects and identifying their types at a larger scale.A-scan, B-scan, and C-scan ultrasonic techniques are commonly used for imaging concrete, with each providing different types of information about the internal structure of the material.A-scan measures the amplitude of the wave received by a single transducer at a specific point, while Bscan provides a cross-sectional image of the material, allowing for the visualization of the location and shape of defects in the plane of the scan [7].C-scan constructs a two-dimensional map of concrete, which is useful for pinpointing the spatial positions of defects [8,9].Other approaches for imaging concrete include coda interferometry using diffusive waves, x-ray computed tomography, and microwave imaging.Planes and Larose [10] showed that simple scattering in concrete occurs near 10-100 kHz frequency, multiple scattering near 100-500 kHz and higher than 500 kHz frequency causes strong absorption.Zhang et al [11] applied the diffuse wave approach to generate 3D imaging of concrete using 16 transmitters and 16 receivers with a frequency range of 80-220 kHz.Zielinska and Rucka [12] developed velocity maps to detect the presence of steel bar or polymer pipe embedded inside the concrete using 18 transmitters and 18 receivers.Achieving the required resolution using linear ultrasonics often requires a significant number of transducers.X-ray computed tomography is a powerful method for assessing the defect evolution inside concrete [10,11]; however, its applicability in the field and large-scale structures is difficult.Microwave imaging requires instrumenting the circumference of concrete with antennas to construct an image using a delay-and-sum beamforming approach [13].The main limitation of microwave imaging is the penetration depth.
The C-scan image is generated by raster-scanning a transducer over the surface of the material and recording the amplitude of the reflected waves at each point.In contrast, ultrasonic imaging provides a three-dimensional view of the material's internal structure, showing a volumetric reconstruction of the material.Ultrasonic imaging is achieved by acquiring multiple C-scans at different angles, typically using an array of transducers, and reconstructing a 3D image of the internal structure from the collected data.The imaging method facilitates mapping internal flaws or material degradation in real-time and in-situ without destructive testing.Many physical waves can be used, for example, visible light, laser, electrical capacitance, resistivity and impedance, x-ray, fluid flow, magnetic induction, and ultrasound [14].
Ultrasonic imaging is most widely used for steel, rock, and concrete.It has several advantages, such as safety, costeffectiveness, efficiency, and rapid application [20].Some examples of ultrasonic imaging of concrete are shown in table 1. Due to the heterogeneity of concrete, the ultrasonic frequency is typically limited up to 100 kHz.In these studies, defects are generally artificially introduced with inclusions during the construction phase.Different properties of ultrasonic waves can be selected to reconstruct the ultrasonic image [21].The most common linear ultrasonic property to generate the image is wave velocity.The resolution of linear ultrasonics is limited by the wavelength of the excitation frequency.
Nonlinear ultrasonic testing (NLUT) is known to be sensitive to defects in sub-wavelength.There are three common NLUT methods: acoustoelastic method, nonlinear elastic wave spectroscopy method, and higher harmonics method.The acoustoelastic method is based on the stress-dependence of the ultrasonic wave velocity [22].The nonlinear elastic wave spectroscopy method is based on the fact that the level of nonlinearity in the elastic response of materials containing Chai et al [19] defects is higher than the one without defects [23].The higher harmonics method is based on the phenomenon wherein the exciting signal is distorted by the nonlinear elastic response of the medium to the incident wave, and the higher harmonic wave is generated [24].The common nonlinear ultrasonic methods applied to the concrete are the second harmonic generation, known as SHG [25], and wave modulation [26].Due to the inherent heterogeneous characteristic of concrete, a typical ultrasonic frequency for the SHG method is limited to below 100 kHz [25].The wave modulation approach tracks the interaction of two monochromatic continuous waves at two separate frequency ranges with defects.The relative strength of two sideband amplitudes is used as an indication of the defects [27].Ciampa et al [28] and Li and Cho [29] demonstrated the tomographic image of steel using nonlinear ultrasonics to detect the presence of corrosion.To the best of our knowledge, no characterization study has been carried out to construct ultrasonic images of concrete based on NLUT.Previous research has primarily focused on either linear or nonlinear ultrasonic methods with little investigation into how to bridge the gap between these two techniques to take advantage of their respective strengths, especially in situations where the size and characteristics of defects are unknown.
In this paper, a measurement framework to generate 2D ultrasonic imaging of concrete using both linear wave velocity and acoustic nonlinearity coefficient obtained from a single measurement is introduced.The framework links the linear and nonlinear ultrasonic parameters to construct 2D ultrasonic images.To demonstrate the dependence of resolution on sensor array and frequency in the linear ultrasonic method, numerical simulations are conducted.The developed framework is evaluated using four calibration samples.The limitations of linear and nonlinear ultrasonic imaging techniques in characterizing concrete defects are described, and their potential solutions to address them are presented.

Linear ultrasonic parameter
The imaging of concrete using linear ultrasonics is based on constructing the wave velocity map.The major measurement variable is the TOF between transmitter and receiver, which is converted into wave velocity using a simple relationship between TOF and propagation length as where L is the propagation distance between transmitter and receiver, V mode is the wave velocity of a particular wave mode.TOF mode is the TOF of wave mode triggered in the structure using an ultrasonic transmitter.For the linear ray-tracing imaging method, TOF is the variable measured, assuming the distances between transmitters and receivers are constant and known.Accurate measurement of TOF mode for detecting the same wave mode for different transmitter-receiver distances is essential to establish an accurate ultrasonic image.In general, the most widely used and simplest arrival time pick-up method is the threshold-based method.The threshold value is pre-defined, and the first-time signal exceeding the threshold is defined as the TOF.The threshold value needs to be set properly to pick up the arrival of the same wave mode.While low computation cost is the advantage of the threshold-based method, it functions less effectively for the signal with a low signal-to-noise ratio and multiple wave modes.To improve the accuracy, several TOF picking algorithms were established in the literature, such as Akaike criteria and Hinkley criteria [30].
An arrival time pick-up algorithm named the outlier-based method was developed recently [31].The outlier is a definition in statistics, which is an observation that lies at an abnormal distance from the other values in a random sample from a population.The outlier-based method to pick TOF is based on the threshold method, which involves selecting an appropriate window and threshold and then calculating the absolute differences in amplitude between neighboring data points within the window.This process creates a new dataset.Typically, the amplitude differences in the new data set conform to a normal distribution at the noise part prior to the onset of the actual signal.The identification of the outliers is predicted on the new dataset's calculated mean and standard deviation, with the first outlier being considered as TOF.The outlier-based method accurately determines arrival time, as shown in two examples in figure 1. Figure 1(a) is derived from a specific instance within the experimental simulation, as elaborated in section 5.In contrast, figure 1(b) shows a sample waveform obtained from the numerical simulation, further detailed in section 4. The arrival times of the experimental and numerical results are approximately 75 µs and 21.8 µs, respectively, due to the size of the numerical model, which is scaled to one-third of the size of the experiment model, resulting in a shorter distance between the transmitter and the receiver in the numerical model.In this study, the outlier method is applied to determine the TOF to calculate wave velocity and wave number.

Nonlinear ultrasonic parameter
The acoustic nonlinearity coefficient using the SHG method, β theory , is proportional to the ratio of the second harmonic amplitude and the square of first harmonic amplitude [32] and calculated as where A 1 and A 2 are the first and second harmonic amplitudes; x is the wave propagation distance; and k indicates the wave number of the first harmonic wave.As the harmonic amplitudes are related to the voltages of harmonic signals extracted from an ultrasonic signal such that β theory is proportional to the voltage amplitudes A 2 v /(A v 1 ) 2 when distance and wave number are constant.The experimental acoustic nonlinearity coefficient is calculated as The most common method to decompose different harmonics is to convert the time domain signal into a frequency domain signal using the fast Fourier transform (FFT) [30,31].However, FFT is not accurate and effective when there is a local material discontinuity in the transient signal [33], potentially resulting in a non-stationary transient signal.The other signal decomposition method is the wavelet transform (WT) [34].Wavelets are mathematical functions that decompose time history data into frequency components.Graps [34] proved that WT could overcome the disadvantages of conventional FFT by analyzing signals containing discontinuities, resulting in greater accuracy and efficiency in signal process applications and defect detection using ultrasonics [8,15,34].
Following the method developed by Mostavi et al [35], the first and second harmonic waveforms are decomposed and employed as inputs for ultrasonic imaging, as illustrated in figures 2 and 3.The parameters extracted from wavelet decomposition to feed into the ultrasonic parameters are demonstrated on a waveform obtained from experimental simulations.Figure 2(a) shows the time history signal and its frequency spectrum.Given the first harmonic transmission signal frequency of 75 kHz, the wavelet spectrogram in figure 2(b) highlights the first and second harmonic frequencies at 75 kHz and 150 kHz.The wavelet spectrogram is a visual representation of a signal's frequency content over time, and it was obtained using complex Morlet wavelet, which is a commonly used wavelet function in wavelet analysis, particularly in the context of time-frequency analysis [36].The wavelet decomposition separates the signal into its frequency components  to calculate the TOFs and amplitudes of the first and second harmonic frequencies.
Before calculating the nonlinearity coefficient, one needs to solve the nonlinear wave equation, which involves higherorder terms that depend on the product of different wave components as shown in equation ( 3).In the process of solving the equation, it is always assumed that the wave number k is a constant for simplicity.However, wave number can vary with different frequencies and wavelengths.In this study, wave velocity extracted from the time domain signal of the first harmonic frequency (f 1 ) is used to calculate wave number using the relationship as ) . ( The schematic of the hybrid measurement approach is illustrated in figure 3. The time domain signals of the first and second harmonic frequencies are extracted from the wavelet spectrogram and plotted in figure 3. The arrival time of the first harmonic signal is determined by the outlier-based method as the TOF input for reconstructing the linear ultrasonic image.The wave number k is obtained using the wave velocity which depends on the TOF of the first harmonic wave as indicated in equation ( 4).The amplitudes of the first and the second harmonic frequencies and k are used to calculate the acoustic nonlinearity coefficient β to construct the nonlinear ultrasonic image.With this hybrid approach, two ultrasonic images of linear wave velocity-based and acoustic nonlinearity-based are generated within a single measurement.The 2D ultrasonic image reconstruction algorithm is presented in the next section.

Linear ultrasonic image reconstruction
The linear ultrasonic image reconstruction includes two methods: ray-based [37] and waveform-based [38].The waveformbased method is widely used in biomedical materials because of its good accuracy and resolution.The ray-based method is typically used for engineering materials due to its advantage of lower computation cost.The outcome of the linear ray-based method is the construction of a velocity distribution map.Straight-line or curved ray methods were previously studied [12], showing good results in the detection of flaws in solids.Herein, the straight ray-based method is applied.The primary experimental input to the algorithm is TOF.The geometry is discretized in a certain pixel depending on the ray paths.The interpretation of pixel is such that the model is divided into a number of uniform grids, where each grid is considered as a single pixel.TOF within each pixel is calculated as where i indicates any ray path; S j represents the slowness within pixel j, and is the inverse of velocity, s m −1 ; n is the total number of rays as a function of the transmitters and receivers; and l ij is the length of ray-path i within a single pixel j.As an example, the ray path in each pixel is illustrated in figure 4.
The layout of transmitters and receivers is defined by green and blue patterns, respectively.The path of wave propagation passing through each pixel is illustrated in figure 4(b).The model is built such that the ray path is considered a straight path between transmitters and receivers, and the slowness, or velocity, in each pixel, remains constant.Multiple wave paths are constructed.The uniform mesh can be created, as shown in figure 4, as N lines × N rows.The pixel index is shown in figure 4(b).To obtain the information of the entire distribution map, equation ( 5) may be extended into matrix form as where m is total number of the pixels, m = N × N in this mesh; n indicates the total number of wave paths depending on the number of the transmitters and receivers.Iterative methods are commonly used to solve the linear equation, which is based on the back-projection updating procedure, algebraic reconstruction technique (ART) and simultaneous iterative reconstruction techniques [39].In this study, the ART method was used to solve equation (5).The algorithm involves two steps: (1) obtaining the [L] matrix depending on the number and placement of sensors, where the TOF measurement serves as input to equation (5); and (2) solving the linear equation using the ART algorithm.The opensource tool (AIR) associated with MATLAB is used [40], where the output variable S, the slowness, inverse of the velocity, provides a velocity visualization within each pixel.The 2D ultrasonic image of the sample is reconstructed accordingly.

Proposed nonlinear ultrasonic image reconstruction
A similar ray-based method is adapted to reconstruct the nonlinear ultrasonic image.The experimental input to the algorithm is the acoustic nonlinearity coefficient in the individual ray path.In acoustic nonlinearity-based ultrasonic imaging, the amplitude of the ultrasonic wave is not linearly proportional to the pressure wave that generates it, due to the nonlinear acoustic properties of the material.Additionally, the heterogeneity of concrete affects the acoustic nonlinearity coefficient depending on transmitter and receiver angle.Therefore, the conventional ART method may not provide sufficient accuracy in ultrasonic imaging based on the acoustic nonlinearity coefficient.Another approach is to use a B-scan technique from two different orthogonal directions, considering only the direct path from the transmitter to the opposite receiver, to reconstruct a surface plot that visualizes the nonlinear ultrasonic image.Section 5 of this paper discuss the specifics of reconstructing the nonlinear ultrasonic imaging and explores the challenges of in-direct way paths.

Description of models
Numerical models using COMSOL Multiphysics software were designed to evaluate the image reconstruction algorithm and understand the resolution of linear ultrasonic imaging.
To reduce the computational time, a 2D plane strain model was selected, representing the cross-section of the cuboid concrete geometry.The excitation signal was a five-cycle sine wave function with a Hanning window as shown in figure 5(a).Three frequencies (75 kHz, 100 kHz and 300 kHz) were modeled to understand the frequency-dependent resolution.
The concrete geometry was a square shape with the length of 75 mm, and modeled with varying size and inclusion forms, see figure 5(b).The length of the inclusion was increased from 15 mm to 25 mm.Examples of received signals are shown in figure 5(c).The TOF was influenced by two factors: the distance between transmitter and receiver, and the presence of an inclusion in the ray path.The transducers were modeled as points and displayed evenly along the side of model.Finite element analysis was conducted using a free triangular mesh with a mesh size of one-twenties of a wavelength.Low-reflecting boundary condition was applied to the four edges.The material properties are summarized in table 2. To evaluate the optimum resolution to detect inclusions depending on the pixel resolution, three models with different transducers were built.The number of pixels in one direction was selected as the same value as the number of transducers because the number of pixels was highly dependent on the selected transducer layout.

Image enhancement
An image post-processing approach through interpolating the points between each pixel was implemented.The cubic interpolation of values at neighboring points in each respective dimension was applied to eliminate the distinct boundary between neighbor pixels.The approach was applied to the image constructed using an excitation frequency of 300 kHz and a 16 × 16 transducer array as an example.As shown in figure 6, image enhancement improves the visualization of inclusions.It is important to note that the method does not affect the resolution and accuracy of the imaging algorithm.Increasing the number of transducers attached to the model can lead to a more precise localization of the inclusion with fewer artificial regions observed outside of its actual position.These findings suggest that increasing the number of transducers can significantly enhance the spatial resolution and   improve the sensitivity to the measurement.The wavelength of 75 kHz frequency in concrete is approximately 40 mm, and therefore, the detectability using 75 kHz is about 20 mm.Consequently, a higher frequency is needed to detect inclusion sizes of 15 mm more clearly while heterogeneous concrete medium limits selecting higher frequency to prevent scattering and absorption.The effect of the proportion of the inclusion size and the base material size was investigated.Although the inclusion may be relatively large, the dimensions of the base material may result in a longer wave propagation depth, which causes more attenuation and changes in the receiver response depending on the transmitter-receiver angle.As shown in figure 8, three sizes of base material as 75 mm, 125 mm, and 225 mm, with an inclusion size of 25 mm, were modeled, resulting in inclusion-to-base material proportions of 1:3, 1:5, and 1:9, respectively.The simulations were conducted using a uniform sensor array of 8 by 8 with an excitation signal frequency of 75 kHz.The results demonstrate that when the base material size is three times that of the inclusion, the presence of inclusion is clearly detected, as evidenced by a significant drop in velocity within the range of the inclusion.However, as the inclusion-to-base material proportion increases, the sensitivity of the method decreases, ultimately resulting in a complete loss of inclusion effect at the inclusion-to-base material proportion of 1:9.This study highlights the importance of  considering the size of the base material in addition to the minimum detectable inclusion size using the linear ultrasonic.

Numerical results
The numerical models were repeated for 100 kHz and 300 kHz excitation frequencies to investigate the effect of the excitation frequency on the imaging results.Figure 9 shows that the detectability of inclusion using 100 kHz and 300 kHz is higher as compared to 75 kHz.Increasing the excitation frequency increases the accurate detection of inclusion with smaller artificial regions.For instance, in the 300 kHz excitation frequency with an inclusion size of 25 mm, a clear indication of the inclusion with minor artificial regions with low wave velocity is observed.In summary, the performance of linear wave velocity-based ultrasonic imaging depends on transducer layout, inclusion-tobase material proportion, and excitation frequency.The number of sensors affects the resolution and the quality of images, although it may be limited by the physically available space to fit the sensors along concrete edges.While higher excitation frequency improves the resolution, heterogeneous characteristics of concrete may cause scattering and absorption of ultrasonic waves and produce inaccurate image reconstruction.

Sample preparation and experimental setup
The ultrasonic image reconstruction algorithm was evaluated using four concrete samples having the dimensions 225 mm × 225 mm × 200 mm, as shown in figure 10.The first sample consisted of pure concrete without inclusions and was used as a benchmark.Expanded polystyrene (EPS) insulation foam with a 25 mm side length was inserted through the concrete section in the second sample.The third and fourth samples had single EPS foam with sizes of 45 mm and 60 mm, respectively.The concrete mixture is given in table 3. The measurements were taken after the specimens were cured in a moisture room for 28 d.
As depicted in figure 10, the measurements were carried out using a PCI-8 data acquisition system and transducers manufactured by MISTRAS Group Inc. R6 (the resonant type with a peak near 60 kHz) and R15 (the resonant type with a peak near 150 kHz) transducers were used as transmitters and receivers, respectively.Transducers were mounted in both x and y directions.Nine receivers and nine transmitters were positioned at 25 mm spacing, with a spacing of 12.5 mm between the first and the nearest edge, as shown in figure 10(a).To mitigate the multiple scattering phenomenon and energy absorption in concrete at the higher wave frequency, the excitation signal was selected as 75 kHz.Moreover, the second harmonic frequency would match the receiver's resonating frequency of 150 kHz, thereby enhancing the amplitude of the second harmonic wave.All transducers were adhered to concrete using hot glue.MATLAB program was used to obtain all the ultrasonic parameters and images.

Linear ultrasonic results
Before building the linear ultrasonic images, the wave velocities along the direct paths of transmitters and receivers' measurements were analyzed to develop B-scan plots.The results of four samples are shown in figure 11.In pure concrete, the wave velocity remained nearly constant with only minor variation, as illustrated in figure 11(a).Conversely, concrete with foam displayed a drop in velocity at the seventh testing point in the x and y direction, as depicted in figure 11(b).However, most paths passing through inclusions showed no  significant changes in wave velocity.Out of the first two concrete samples, it is apparent that using a limited amount of transducer and linear ultrasonic with a 1:9 proportion of inclusion-to-base material proportion did not perform well in detecting defects, similar to the results obtained from the numerical models.However, when the size of the foam inclusion was increased to 45 mm and 60 mm, which are equivalent to 1:5 and 1:3.75, there is a clear drop in the wave velocity in the path passing through the inclusion.A consistent trend is observed in both the x and y directions.One important point to note is that the use of slag in the concrete mixture and a lower water-to-cement ratio resulted in higher concrete strength f ′ c in the range of 45 MPa to 48 MPa (6500-7000 psi) for these concrete samples.The corresponding Young's modulus is in the range of 32 000 MPa to 33 000 MPa (4600 ksi to 4800 ksi).As a result, the wave propagation velocity in this high-performance concrete sample is relatively higher as compared to the value used in numerical models (Young's modulus taken as 25 000 MPa).This is a significant indicator that ultrasonic imaging should focus on the variation of relative wave velocity as it depends on concrete strength.
To expand from a 1D B-scan plot to a 2D linear ultrasonic image, the wave velocities through the angled ray paths were extracted.Ultrasonic images were constructed, shown in figure 12, using the method described above.Due to the limited number of sensors, excitation frequency, and low inclusion-tobase material proportion, the linear wave velocity-based image could not accurately detect the position of the small EPS foams as expected.The artificial regions present in the pure concrete sample can be attributed to the limitations of the ray tracing algorithm.In figure 12(b), low wave velocity regions coincide with the positions of the inclusion: from the left bottom corner to the right upper corner.However, the positions of the inclusion cannot be discerned if they were not known beforehand.
As depicted in figures 12(c) and (d), the linear ultrasonic image of the larger inclusion can be detected well due to higher inclusion-to-base material proportion for the selected frequency and transducer array.The inclusion in the left corner of the concrete is clearly visible, with a significant drop in wave velocity, consistent with the result obtained from the numerical simulation.The inclusion-to-base material proportion is an important factor impacting the sensitivity of linear ultrasonic testing, as evidenced by these findings.

Nonlinear ultrasonic results
As compared to the linear ultrasonic results, the B-scan results of the nonlinearity parameter have distinct changes at the coordinates of inclusions.As shown in figure 13(a), the nonlinearity parameter appears with a minor variance in the pure concrete.The value of the nonlinearity parameter rises when the wave propagates through foams in both orientations.The nonlinear ultrasonic results show better performance in detecting inclusions as compared to the linear ultrasonics.
There are several factors that could potentially influence the value of the nonlinearity coefficient.Upon comparing the samples with varying sizes of inclusions, it becomes apparent that increasing the size of the inclusion leads to greater material discontinuities, resulting in a higher acoustic nonlinearity coefficient.In addition, the inclusion-to-receiver distance is found to be another influential factor.As the inclusion-to-receiver distance increases, the nonlinearity coefficient decreases.This phenomenon can be attributed to the fact that the acoustic nonlinearity coefficient is proportional to the ratio of the second harmonic amplitude to the square of the first harmonic amplitude.After the generation of the second harmonic signal via wave distortion at the inclusion, it propagates toward the receiver with higher energy absorption and attenuation than the first harmonic signal due to the multiple scattering phenomena of high frequency signal in concrete.On the other hand, the first harmonic experiences a lower attenuation throughout the wave propagation distance.Therefore, the longer the distance that second harmonic travels, the lower acoustic nonlinearity that we observe.
Consistent trends are observed across all the wave paths with the inclusion, as shown in figure 13.By utilizing this trend, it is possible to estimate the relative location of the inclusion by comparing it with other inclusions or by estimating the position of a single inclusion.As shown in figure 13(b), the inclusion located closer to the receiver exhibits higher acoustic nonlinearity in the x-direction, as represented by x 7 .Similarly, the inclusion located farther away from the receiver along the y-direction is denoted by a lower value of y 7 .This capability enables the extension of the 1D scan to 2D by estimating the position of the location using the data obtained from the other orientation.
Expanding the 1D B-scan image to a 2D ultrasonic image using the acoustic nonlinearity coefficients of angled ray paths raises some challenges.The nonlinear ultrasonic coefficient in concrete can be direction-dependent, which means that the wave behavior varies depending on the angle between the wave front to the defect and the distance between the transmitters and receivers.For instance, if the presence of the defect remains in the same position, but the ultrasonic wavefront encounters the defect at a different angle, the waveform may vary [32].When the wavefront is perpendicular to the defect, the interaction between the waves and the defect is expected to be stronger, which leads to higher acoustic nonlinearity coefficients.Additionally, the distance between the transmitter and receiver can also impact the nonlinear ultrasonic measurements, as discussed above.To explain the direction-dependent characteristics, the nonlinearity coefficients variation of one transmitter to the other nine receivers of the pure concrete sample and the sample with two foam inclusions are analyzed, as shown in figure 14.
When comparing each path in pure concrete without any inclusions, it is noticeable that the nonlinearity coefficient varies depending on the angle of the waves, which results in varying distances between the transmitter and receiver.With the exception of the first two paths at the bottom, a decreasing trend is observed as the travel distance increases.Figure 14(b) shows that when a wave passes through a single foam inclusion, such as path 5 and path 6, there is an increase in the nonlinearity coefficient.In paths 8 and 9, a significant decrease in the nonlinearity coefficient is observed when the wave passes through two inclusions.This is attributed to the absorption of the second harmonic wave generated within the first inclusion by the second inclusion, leading to a reduction in energy.The study found that this trend was consistent across all analyzed ray paths.Therefore, selecting appropriate ray paths is crucial for constructing accurate nonlinear ultrasonic images that account for the effect of multiple inclusions on the energy   absorption of the second harmonic wave.Current ray-tracing algorithms do not account for the direction-dependent characteristics and multiple-inclusion behavior, which can lead to inaccurate results if all the ray paths are included in the algorithm.To overcome these limitations, the use of the 2D surface plot combining x and y orientation B-scan image is proposed to construct a more reliable and effective nonlinear ultrasonic image.
The 2D ultrasonic image reconstructed using the acoustic nonlinearity coefficient is depicted in figure 15.The nonlinearity coefficient varies significantly with the inclusion size.Therefore, the scale of the image has been adjusted to improve visualization.Pure concrete exhibits a minor acoustic nonlinearity coefficient, as shown in figure 15(a).In contrast, when two small EPS foams are present in the concrete, a clearer image that differentiates the regions of inclusions is obtained.The positions of two foams can be clearly detected, along with two false detections.These false detections are caused by multiplying the B-scan image in both the x and y directions, resulting in a total of four peaks.As discussed in section 5.3, the nonlinearity coefficients are influenced by the distance between the inclusion and the receiver.
If multiple inclusions are present within the concrete, the nonlinearity coefficients of acoustic signals can be used to determine their true or false positions.For example, in figure 13(b), the nonlinearity coefficient of x 7 is higher than that of x 2 , indicating that the inclusion along x 7 is located closer to the right edge of the concrete (the position of receivers) than the inclusion along x 2 .Similarly, the nonlinearity coefficient of y 2 is higher than that of y 7 , indicating that the inclusion along y 2 is located closer to the bottom edge (the position of receivers) of the concrete than the inclusion along y 7 .By comparing the nonlinearity coefficient of two peaks in the x and y direction, the inclusion-to-receiver distance can be estimated, allowing for the identification of false positions when the 2D image is constructed using the direct path results, as shown in figure 15(b).Figures 15(c) and (d) show that larger foam inclusions in a single position can be clearly identified using a limited number of transducer array.The region with brighter color indicates a higher nonlinearity coefficient, proving that larger inclusion exhibits a stronger nonlinear response.These results demonstrate that the nonlinear ultrasonic method can detect low-density inclusions and localize them accurately with high resolution, even when a limited number of transducer array is used.

Conclusions
An integrated measurement of linear and nonlinear ultrasonic imaging was developed to detect sub-surface defects in concrete structures.The same sensor array was used for both measurements while the received signals were decomposed into harmonics to generate inputs for linear wave velocitybased imaging and acoustic nonlinearity coefficient-based imaging.The image reconstruction models were developed under the assumption that the ray path was straight between transmitters and receivers.The linear wave velocity image was developed such that total slowness in each path was the summation of the slowness within each pixel region that the path passed.It was numerically shown that the image resolution of linear wave velocity-based images was dependent on the number of transmitters and receivers, which was limited by the physical available space around each concrete specimen and the size of the sensors.Increasing the transducer array, such as a 16 × 16 configuration with 5 mm sensor spacing, enhances image resolution and the detectability of defects.The sensitivity of linear ultrasonic imaging was also affected by the inclusion-to-base material proportion.For a sensor array of 9 × 9 with 25 mm sensor spacing, the linear ultrasonic image was effective when the inclusionto-base material ratio was high (e.g.1:3), but its detection capability decreased significantly when the ratio was lower (e.g.1:9).The acoustic nonlinearity coefficient image was reconstructed using the surface plot, which were calculated by multiplying the nonlinearity coefficient in x and y orientations as the conventional ray-tracing method does not account for the direction-dependent characteristic nature of concrete.Combining linear and nonlinear ultrasonic images enabled the identification of defects at multiple scales, leveraging the strengths of each method.The nonlinear ultrasonic image could detect smaller defects than the linear ultrasonic image when the same number of sensors was used.For the same sensor array of 9 × 9, the nonlinear ultrasonic image could detect the inclusion-to-base material ratio down to 1:9.However, their high sensitivity to material discontinuities and experimental variability poses drawbacks.In contrast, linear ultrasonic can achieve higher image resolution, taking advantage of using all the ray paths included in the image reconstruction algorithm.Combining both linear and nonlinear ultrasonic imaging within a single measurement provides a more reliable assessment, especially when the defect size and position are unknown.Future work will involve developing reconstruction algorithms that account for the direction-dependent characteristics and designing a smaller sensor array to increase the spatial resolution of ultrasonic images.Moreover, the numerical simulation of nonlinearity in concrete structures will be studied.

Figure 1 .
Figure 1.The arrival times of two waveforms highlighted with a red line by the outlier-based method, (a) the arrival time of the experimental waveform and (b) the arrival time of the numerical waveform.

Figure 2 .
Figure 2. (a) The received signal in the time domain and its frequency spectrum, (b) the wavelet spectrogram indicating 1st harmonic and 2nd harmonic frequencies.

Figure 3 .
Figure 3.The schematic for calculating inputs for linear and nonlinear ultrasonic imaging, TOF for linear ultrasonic imaging and β for nonlinear ultrasonic imaging.

Figure 4 .
Figure 4. (a) Schematic of straight ray path on a two-dimensional geometry and (b) the regulation of pixel index.

Figure 7 (
Figure 7(a) shows three models with varying inclusion sizes.Figures 7(b)-(d) show the linear ultrasonic images of transducer layouts for 8 × 8, 12 × 12, and 16 × 16, respectively, with an excitation signal frequency of 75 kHz.The red dashed box highlights the inclusion location while the other regions with low wave velocity is artifacts of the algorithm.The region of low wave velocity observed in the 25 mm inclusion is more focused than that in the 15 mm inclusion as detecting larger inclusions with the linear ultrasonics becomes easier.The 12 × 12 and 16 × 16 transducer arrays clearly detect inclusion sizes of 20 mm and 25 mm.In contrast, the sample with 15 mm inclusion results in a fuzzy image with numerous artificial regions with low wave velocity that could impede the accurate identification of the real inclusion.Increasing the number of transducers improves the resolution of imaging, even when the inclusion size is the same, as depicted in figures 7(b)-(d).Increasing the number of transducers attached to the model can lead to a more precise localization of the inclusion with fewer artificial regions observed outside of its actual position.These findings suggest that increasing the number of transducers can significantly enhance the spatial resolution and

Figure 5 .
Figure 5.The numerical model details (a) the excitation signal in time and frequency domains, (b) 2D concrete model with inclusions indicating transmitters T 1 to Tn and receivers R 1 to Rn and (c) examples of received waveforms along different rays.

Figure 6 .
Figure 6.An example of linear wave velocity-based imaging of concrete with foam inclusion highlighted with red dashed line, (a) original image and (b) enhanced image.

Figure 8 .
Figure 8.The simulation results of 75 kHz using an 8 × 8 transducer array (a) the dimensions of inclusion and base material, (b) the proportion of 1:3, (c) the proportion of 1:5, and (d) the proportion of 1:9.

Figure 9 .
Figure 9.The influence of excitation signal on the imaging result using a 16 × 16 transducer array and three inclusion sizes, (a) the profiles of models, (b) the ultrasonic images using 300 kHz excitation signal, (c) the ultrasonic images using 100 kHz excitation and (d) the ultrasonic images using 75 kHz excitation.The black circles indicate the artificial regions with low wave velocity.

Figure 10 .
Figure 10.Concrete samples, (a) no inclusion, (b) two EPS foams of 25 mm, (c) EPS foam of 45 mm, (d) EPS foam of 60 mm is instrumented using transmitter and receiver array, which are connected to the (e) data acquisition system.

Figure 11 .
Figure 11.The wave velocities obtained from the B scan of transmitter-receiver pairs, (a) pure concrete, (b) concrete with two EPS foams, (c) concrete with EPS foam of 45 mm, and (d) concrete with EPS foam of 60 mm.

Figure 12 .
Figure 12.The experimental linear velocity-based ultrasonic images of four concrete samples, (a) pure concrete, (b) concrete with two EPS foams (c) concrete with a foam inclusion of 45 mm, and (d) concrete with a foam inclusion of 60 mm.The red dashed boxes highlight the actual positions of inclusions.

Figure 13 .
Figure 13.The acoustic nonlinearity coefficient obtained from the direct paths of transmitters (T) and receivers (R) along both x and y directions, (a) pure concrete, (b) concrete with two EPS foams, x 2 , x 7, y 2 and y 7 represent the acoustic nonlinearity in the second and seventh paths in x and y directions, respectively (c) concrete with EPS foam of 45 mm and (d) concrete with EPS foam of 60 mm.

Figure 14 .
Figure 14.The acoustic nonlinearity coefficient obtained from one transmitter to other nine receivers in (a) pure concrete sample and (b) concrete with small EPS foams.

Figure 15 .
Figure 15.The 2D ultrasonic image based on the acoustic nonlinearity coefficient (a) pure concrete, (b) concrete with two EPS foams, (c) concrete with EPS foam of 45 mm, and (d) concrete with EPS foam of 60 mm.The red dashed boxes highlight the actual inclusion positions.

Table 1 .
Applications of ultrasonic imaging to concrete.

Table 2 .
Material properties used in numerical simulations.