Abstract
Radar detection has proven to be an effective, nondestructive test for the determination of the quality of wood-based materials, especially in the wooden structures of ancient buildings and trees. However, the results are usually inaccurate, and it is difficult to interpret internal anomalies due to the moisture content of wood, individual differences, and other factors. In this paper, a new measurement method is proposed based on the use of ground-penetrating radar (GPR) for abnormality localization and imaging. Firstly, the time delay of the reflected signal in the inner trees is analyzed with matched filter and Hilbert detections. Secondly, the two approaches are compared with the use of a forward model, and the Hilbert algorithm is found to be more accurate. Thirdly, a laser scanner is used to collect contour data and determine the location and characteristics of internal tree anomalies. Lastly, the proposed method is tested on ancient willows at the Summer Palace. The results show that the error in the depth and area estimates of the anomalies was within 10% and 5%, respectively. Consequently, the GPR method for locating the anomalies in trees is feasible, and a laser scanner combined with contour data can present the size of the abnormal regions within the trees.
1. Introduction
The preservation and maintenance of ancient trees are global concerns that are related to the protection of ancient cultures. Many wild trees are endangered by natural defects and biological degradation. In urban, trees around roads are easily broken and collapse in adverse weather conditions such as strong winds and heavy snow storms, bringing unnecessary losses to lives and property. Therefore, it is very important to have nondestructive testing methods for wood to assure that trees grow well and that the famous ancient trees and buildings are well preserved.
In years past, the evaluation of wood was mainly by visual inspection. Additional methods have gradually developed, and methods currently include stress wave [1], acoustic pulse [2], tree needle [3], ultrasonic [4], X-ray scanning [5], radar, and other methods. Most methods have some disadvantages, such as being time-consuming, invasive, and complex. Compared to other nondestructive detection technologies for wood, radar detection with radio frequency waves is totally nondestructive, convenient to operate, highly efficient, and environmentally friendly, making it a practical method [6].
The radar is a well-established nondestructive testing (NDT) technique that uses radar waves to detect objects and determine their distance based on the echoes they reflect [7]. Unlike traditional radar systems, ground-penetrating radar (GPR) systems are mainly used to detect and measure the depth of inhomogeneity (either defects or layers) in a dielectric medium. Detection is achieved by comparing the power of the scattered electromagnetic (EM) waves produced by the contrast in the dielectric properties of the medium and the inhomogeneity to a preset threshold above the receiver noise level [8]. GPR is a NDT method for locating and accessing structural objects, including pavement material layer thicknesses and properties [7, 8]. This method is currently used in many areas as an investigation tool in the protection of cultural heritage [9], pavement detection [10, 11], soil layers [12], and civil engineering [13, 14].
In recent years, the acceptance and use of GPR tree analyses have gradually grown. Several types of measurements have been made to analyze the use of GPR to detect tree root systems, with which it was possible to detect medium and large roots with groups of small roots. In the study, two antennae, 900 MHz and 1 GHz, were used to detect tree root systems by applying GPR polarimetry techniques [15]. Trunk measurements were made using GPR to sense vertically along the trunk with a metal shield for estimation of the average of tree trunks. In the reflection mode, the transmitter and receiver are within the same box [15]. Other studies have also used radar waves to study the dielectric properties of trees [16–19]. Fu et al. reported a ray-based tomography method using GPR to reconstruct the internal structure of the trunk cross-section of a living oak tree. The point cloud technique is used to depict the shape and size of the trunk cross-section [20]. Butnor et al. used GPR to estimate decay volumes in living coniferous trees, and the threshold method was used to detect the defects. The defected-area error was approximately 18.9% [21]. Halabe et al. showed that GPR can be used to accurately identify subsurface defects such as knots, decay, and metallic nails inside logs that were not visible from outside observations. Combining the use of dielectric properties and the threshold method, the depth error detected at the sawmill was approximately 10% [22]. Lin et al. reported a 17.78–52.61% error in the tree trunk decay area using a two-dimensional stress wave method for detection [23]. To develop an automatic detection system for wood defects, Devaru and Gopalakrishnan studied the algorithm for detecting defects in logs using GPR [24]. Perez-Gracia et al. used GPR in buildings to locate damaged beams and identify internal structures. The study found that differences in reflections due to damaged beams were clearly highlighted in GPR images, while differences were less apparent in healthy beams [25]. Overall, GPR has proved to be a useful tool for detecting tree trunk and wood abnormalities.
Although GPR technology has developed for detection of tree trunk abnormalities, its routine application for evaluation of wood remains minimal and incomprehensive. The limited use of GPR for tree evaluation is due mainly to the lack of automated procedures for data analysis, relatively imprecise accuracy in defect location, and difficulties illustrating the GPR data collected during the tree analysis.
In this paper, both GPR and laser scanning radar are applied to accurately locate internal abnormalities in trees and present the size of the abnormal regions using four steps. First, the desired reflections of the GPR signals and contour information are collected; second, the reflection time delays are determined to develop accurate estimates; third, the medium’s dielectric properties, which are used to estimate the propagation speed of EM waves within the tree, are gauged; and finally, the contour and location information are combined to present the size of the abnormal regions within the trees.
2. Materials and Methods
2.1. Radar Data Acquisition Methodology
This research was conducted using the TRU (Tree Radar) System from Treewin Corporation in America. As shown in Figure 1, the system consists of two parts, a SIR 3000 GPR unit equipped with a 900 MHz center frequency antenna from GSSI (Geophysical Survey Systems Inc., USA).
Figure 2 shows the procedure of the radar data acquisition during the study of Salix babylonica L at the Summer Palace in Beijing of China.
For collecting data, the antenna was placed against the trunk of a tree and moved circumferentially to acquire a complete radar gram of the selected trunk elevation. Every 5 mm, a reflection trace or waveform comprised of the amplitude of reflected energy and the two-way trace time was collected. The traces which contain 255 samples were stacked to create a radar gram.
2.2. Contour Information Acquisition
The structure of the trees in this study was mostly irregular columnar. The cylinder and irregular cylinder approximation methods were applied to obtain the outer contour of the trees to achieve the accurate location of the internal condition of the trees. The cylinder cross-sections were transformed to polar coordinates, making them approximately circular for use in the analysis with software. For irregular columnar-shaped trees, three-dimensional coordinates were established based on the three-dimensional laser scan of the outer contour. The three-dimensional laser scanner contains a laser-emitting device and a laser-detecting device. Two stepper motors control the laser beam movement in horizontal and vertical directions. The scanner emits a laser beam that reflects off an obstacle. When the laser detector in the scanner detects reflected light, measurement of the scanned spot is complete [26, 27]. By measuring the time difference between transmission and reception of the laser, the laser propagation distance can be calculated. Thus, three-dimensional coordinates for objects can be determined to describe the surface of the target, such as the outside contour of trees [28, 29].
2.3. Theoretical Background
As explained above, radar uses radio waves to detect objects and determine their distance based on the echoes they reflect [7]. The GPR is a special type of radar that detects different electrical properties within the medium it is penetrating [30]. Several studies have illustrated the capability of GPR to study wood and detect internal reflections caused by humidity, damage, or disease. Using the GPR, a transmitting antenna emits short pulses of the electromagnetic wave into the trees. These pulses are partly reflected when they encounter an anomaly (i.e., a void or decay) with different dielectric properties and partly transmitted into deeper layers. The arrival time and amplitude of the reflected wave vary with the location and the dielectric discontinuities of material (e.g., air/bark and bark/anomaly) [21]. Furthermore, the difference of the dielectric properties of the medium can also affect the resolution and accuracy [31].
Prior studies have established the use of a forward model of radar electromagnetic wave propagation for layered structures to verify the expected propagation of the reflected radio waves. The forward model is composed of three layers of materials, the surface layer (ε = 5, cm), the intermediate layer (ε = 13, cm), and the bottom layer (ε = 3, cm). Figure 3 presents the characteristics of the reflected signal. The Riker wavelet was chosen as the incident pulse in the forward model (900 MHz, Figure 4(a)).
(a) The incident pulse is Riker
(b) The echo wave of the forward model
In such a case, the electromagnetic wave is shot in the vertical direction. According to the Fresnel formula [3]: (1) the difference in the electromagnetic properties of the medium on both sides of the interface and the intensity of the reflected wave are positively correlated and the properties of the medium on both sides can be deduced by the amplitude of the reflected waves; and (2) when waves pass from a low to a high dielectric constant medium, the reflection coefficient is negative, indicating that the direction of the amplitude of the reflected wave is the reverse of the incident wave. In contrast, when waves pass from a high to a low dielectric constant medium, the reflection coefficient is positive, indicating that the direction of the amplitude of the reflected wave is the same as the incident wave. For example, the amplitude direction of the reflected radar wave is the reverse of the incident wave at the interface of the air and the trunk. In contrast, the amplitude direction of the reflected radar waves at the interface of a tree cavity and the trunk is the same as the incident waves. The reflections from tree bark and the reflections from hollow areas are in the opposite directions. Therefore, the amplitude and direction of reflected radio waves are important bases for differentiating interfacial electronic properties. The inverse of the echo waveform between the 1/2 layer and 2/3 layer is shown in Figure 5(b). It has verified the judgment basis.
(a) The flow chart of the matched filter method
(b) The detected layers of the single-channel wave using matched filter
To accurately locate abnormalities within a tree, the time delay of the reflections needs to be calculated accurately, and the medium’s dielectric constant is used to estimate the propagation speed of the EM waves within the tree.
Consequently, the GPR signal in the time domain is assumed to be composed of a series of scaled and time-delayed replicas of the incident pulse , as indicated by the following [32]: where is the incident GPR pulse, is the total number of layers composing the pavement system and is the relative amplitude of the reflected pulse at the th interface, and is additive noise and is the two-way travel time through the th layer. The intensity of the recorded reflections is proportional to the strength of the contrast in the dielectric properties between the layers. The amplitudes of these reflections, called and , are used for the computation of the layer properties [8]: where ε is the dielectric value of the surfacing layer, is the amplitude of the reflection from the surface, is the amplitude of the reflection from the intermediate layer, and Am is the amplitude of the reflection from a large metal plate. Since metal is a good conductor, it can be considered a perfect EM reflector. can be the amplitude of the incidental GPR signal, which is determined by collecting GPR data over a large and flat metal plate placed on the pavement surface [7].
Echo delay time can be used to determine the depth of the target. The two-way travel time, , from emission to detection is measured, and distance of the target from source, , can be determined using the following equation: where is the wave velocity in the medium, is the dielectric constant of the medium, is the velocity of electromagnetic waves in a vacuum, and is the depth of the target of interest.
3. Proposed Methods of Time-Delay Estimation
3.1. Matched Filter Detector
A matched filter allows a maximum signal-to-noise ratio (SNR) [8]. The impulse response of the matched filter is presented as follows: where is the duration of the signal . During the time , the matched filter output exhibits a maximum SNR. A signal is declared to be present with the method when the corresponding matched filter output exceeds the fixed threshold, , at time T. Under these conditions, the matched filter output has a maximum (or a minimum, depending on whether the polarity of the detected signal is reversed compared to the original signal) at time . The threshold, , used with the matched filter is determined from the maximum tolerable probability of a false alarm, , based on the Neyman-Pearson criterion [8, 33]:
The matched filter algorithm for detecting multiple reflected pulses is summarized in Figure 5(a).
As shown in Figure 5(b), we adopted the matched filter method to analyze the signal channel echo data using a forward model. The four inset pictures from top to bottom are, respectively, the single-channel echo signal and the reflected signals of the first, second, and third layers extracted from the echo signal.
3.2. Hilbert Detector
The traditional Fourier transform is suitable for smoothing linear signals or a periodic signal and converting the signal to obtain spectral information on the overall situation. However, in engineering applications of radar, the collected echo data is typically a nonstationary and nonlinear signal. For nonstationary and nonlinear signals, time-frequency analysis is usually used to obtain both the short-time Fourier transform and wavelet transform of the signals. The basic functions of the time-frequency analysis method based on the Fourier transform are relatively fixed. The analysis of wavelet transformation requires selection of the appropriate basic functions to achieve the best results. Thus, these time-frequency analysis methods are not very adaptable. When these methods are used to process a signal, they are prone to showing unwanted signals. Huang et al. proposed a new method of time-frequency analysis using the Hilbert-Huang transform [34]. This new method for the analysis of a nonstationary, nonlinear signal is more intuitive and highly adaptive. Using this approach, the signal is processed by empirical mode decomposition (EMD), from which we obtain a finite number of intrinsic mode components (IMF) and a residual signal that represent the trend in signal change. Each IMF obtained by using Hilbert transform time is put forward for frequency analysis [35].
According to the Rosenfeld subband product theory, if the noise signal can be decomposed into multiple frequency subband components, each component of the point can highlight the useful signal, and the noise can be eliminated by point product operation. Based on this view, K. H. Kim and S. J. Kim [36] proposed using a wavelet detector to detect biomedical signals at a low SNR. Inspired by the wavelet detector, in this paper, IMF components are used to construct a detection algorithm. The Hilbert algorithm [37] includes the following four steps (Figure 6(a)): (1)Process the echo data by EMD and obtain the corresponding IMF (2)Calculate the dot product of the absolute value of all components IMF (3)Because the echo signal contains the positive and negative peaks, the received signal is composed of the multipeak pulses, easily leading to false detections. The signal need to be smoothed by window function. Convolute using and window functions to complete smoothing (4)Compare to the threshold . If exceeds the threshold, take the time delay of the maximum value as the time delay of the detected pulse. Otherwise, there is no detected targets
(a) The flow chart of Hilbert algorithm
(b) The detected layers of the single-channel wave using Hilbert algorithm
As shown in Figure 6(b), we adopted the Hilbert filter method to analyze the signal channel echo data using a forward model. The reflected signals of the first, second, and third layers extracted from the echo signal are presented. The Hilbert product algorithm has been verified to be able to locate different media layers.
3.3. Performance Comparison of the Detectors
We use a forward model test based on GPRMAX to test the accuracy of the two algorithms [38]. GPRMAX is based on FDTD (finite-difference time-domain) and PML (perfectly matched layer) algorithms to build a model. A wedge forward model was used to compare the above algorithms, as shown in Figure 7(a). The model is 50 cm in length (: 5 cm) and 110 cm in width (:10 cm). In this process, 36 different points are selected, and the error of each point is analyzed by a Monte Carlo method radar scanning on the forward model from right to left and top to bottom, yielding 36 sets of experimental data. An error analysis was conducted on each set of data. The distance between the surface and the intermediate layer was estimated for each of the two different detection algorithms. Moreover, the accuracy of the surface and the intermediate layer thickness estimates were known based on the wedge forward model. The accuracy of the two different algorithms for layer recognition was calculated with the following equation:
(a) Forward model of wedge
(b) Result of comparison
As seen in Figure 7(b), the Hilbert detector has a higher resolution than that of the matched filter. The A-line in the graph represents the interface between the first layer and the second layer. The B-line in the figure represents the computed interface that can be distinguished by a matched filter. The C-line in the figure represents the computed interface that can be distinguished by the Hilbert detector. Through data calculations, the minimum distance to the recognized layer using the matched filter and the Hilbert detectors was determined to be approximately 7.3 cm and 6.7 cm, respectively. Thus, the Hilbert detector has a higher recognition resolution.
4. Experiment and Data Analysis
4.1. Analysis of Field Data
The GPR technique makes use of electromagnetic waves to recognize decay in trees. The subsurface interfaces of materials with different dielectric constants can be identified based on the amplitudes of reflected waves. The scanned data containing these reflections can be analyzed to differentiate defective and normal areas of logs. These pulses are reflected from subsurface interfaces or boundaries between materials with different dielectric constants (interface between good and defective wood) and received by the antenna. This study adopted the Hilbert detector to estimate the time delay of the echo data to more accurately locate the position of the target. The time delay of the echo data can be used to deduce the location of the anomalies in the wood. For irregular columnar-shaped trees, three-dimensional coordinates were acquired using a three-dimensional laser scanner to track the outer contours. The location of the anomaly was combined with the contour data to determine the size of an anomaly within a tree.
This experiment was conducted on Salix babylonica L at the Summer Palace. Image analysis was conducted on three samples from two trees, T1 and T2 (Figures 8(b) and 9(b)). After a series of radar data preprocessing steps, the Hilbert detection algorithm was used to locate inner anomalies (decay or cavities) in the trees. The raw data was preprocessed to amplify the target signal, remove bad traces, and replace them with computed from their near neighbors. Figures 8(a) and 9(a) show B-scans around the samples using the 900 MHz antenna. The results of this analysis indicated that the area enclosed in the dashed green line was a possible location of air/bark. The dotted red line indicated more certain anomalies (voids, decay). As shown in Figures 8(d) and 9(d), combining the location of anomalies with contour data identifies the size of the abnormal regions within the trees. The cross-sections from the cut tree are shown in Figures 8(b) and 9(b). There is a strong correlation between changes in the GPR-scan signal characteristics and the actual location of anomalies observed in the sawed samples.
(a) B-scan
(b) T1
(c) Analysis of Treewin
(d) Abnormal dimensional analysis chart
(a) B-scan
(b) T2
(c) Analysis of Treewin
(d) Abnormal dimensional analysis chart
In the process of experimental data analysis, the calculated data about the depth distance were all derived from formulas (4) and (5) in Section 2.3. As for T1, the distance from the bark to the cavity was approximately 4 cm, while the software analysis estimated this distance as 3.91 cm. Compared to the tested data, the error rate between the real and estimated values was approximately 2.25%. The distance of abnormal point A, B, C, D, and E from the bark in the B-scan figure was calculated by (3). For T2, the actual distance of the anomalies from the bark was measured by drilling into the cross-section samples (A, B, C, D, and E). We then compared the software-estimated thicknesses to physical measurements from the cores. The calculated depth error and error variance are shown in Table 1.
The actual abnormal area of samples can be calculated using a grid method. The same area was estimated using the Treewin Corporation threshold to analyze the radar data combined with the polar coordinate conversion method to develop a two-dimensional image of the inner anomalies in the trees. The estimated areas using Treewin software analysis had an error rate of approximately 20% compared to the actual (physical) measurements with the grid method. A third approach, based on the location of the anomalies in the trees using Hilbert detection and combined with tree contour data obtained by the three-dimensional laser, introduced contour tracing and a coordinate transformation method to generate two-dimensional images of anomalies within the trees. The results show that the estimated area values of that method have an error rate of approximately 5% compared to that of the actual measurements from the grid method. Table 1 shows the area error ratio, S2 for T2. The grid method calculates the area by dividing the irregular pattern into a number of squares with an area of 1 and plotting the area with an edge area of less than one. Finally, the area of the irregular figure is obtained by calculating the number of squares.
In summary, it can be seen from Table 1 that the use of Hilbert algorithm with 3D scan data and the contour tracing method is able to reduce the error in tree abnormality depth estimates by approximately 10%. In the two-dimensional images shown in the hierarchical diagram, the error in abnormality areas calculated with the Hilbert approach is approximately 5% compared to that of the physical measurement-based cross-sectional areas. This result is an improvement in the accuracy of abnormality area estimates compared to other NDT approaches.
4.2. Point Cloud Technique
In the paper, we used commercial software (Agisoft Photoscan) to acquire point cloud data [20]. A point cloud is a set of data points in a three-dimensional coordinate system, usually defined by X, Y, and Z coordinates. Digital photographs of the trees are the input for this software. The input pictures do not require any preprocessing by the user. Figure 10(a) shows the 3D simulation diagram, Figure 10(b) shows a photograph of the tree, and Figure 10(c) shows cross-sections of the tree at the three heights: 0.9 m, 1.2 m, and 1.5 m. The black area of Figure 10(c) is the internal structure of the cross-section of the tree.
5. Conclusions
This paper proposed a method that combined GPR and laser scanning contour data for locating internal abnormalities in trees. The method determines the size of the abnormal regions. Matched filter and Hilbert detector approaches were applied to obtain the time delay of radar echo data, amplitude ratios were used to deduce the dielectric constant of medium layers, and the forward model (FDTD: finite-difference time-domain) was introduced to verify the algorithm. Finally, the Hilbert algorithm was found to have higher accuracy than that of the matched filter algorithm. Using the above methods, samples of Salix babylonica L at the Summer Palace were analyzed, and the results showed that the error in abnormality depth and area can be reduced to within 10% and 5%, respectively.
This approach greatly improves the accuracy of tree abnormality measurements compared to the Treewin Corporation analysis system. The innovations are as follows. (1) The Treewin Corporation method detected the abnormality position within the tree, using the threshold method to develop a two-dimensional image of anomalies in the trees. In this paper, the abnormal area estimates obtained by the Hilbert method were more accurate than the estimates based on the matched filter method. Hilbert detection, combined with the coordinate transformation method, was used to develop a two-dimensional visualization of anomalies within trees; combined with contour information from laser scanning, the visualization of irregular-shaped trees was achieved, and the accuracy of the abnormality area estimates within the trees was improved. (2) The Treewin system analysis of samples uses a fixed dielectric constant of 13 [21], but the actual dielectric constant of each layer in the tree was slightly different. In this paper, the relative dielectric permittivity of the trunk was obtained from the amplitude of the compared waves, as shown in (2), which resulted in more accurate estimation of the internal abnormal area. (3) In the Lorenzo et al. study of tree detection with GPR, the detection of the trunk was vertical, and the antenna moved along the trunk, so the cross-sectional profile of the trees could not be obtained [15]. In this paper, a 360-degree scan of the trunk is made at selected heights, providing information on the cross-section of the trunk and allowing the generation of two-dimensional cross-sectional visualizations. In the future, this approach can also be used for transverse measurements at multiple levels, allowing generation of three-dimensional stereoscopic images of trees and reflection of internal tree information in a more intuitive manner within certain heights.
Presently, radar cannot effectively detect thin layers due to poor resolution in the time domain. On the premise that the hardware limits this detection, methods can be introduced such as MUSIC (multiple signal classification) [29] and other high-resolution frequency-domain approaches that can detect wooden bodies, allowing further analysis of thin layers of wood.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research is supported by the National Natural Science Foundations of China (Grant No. 31600589) and the Fundamental Research Funds for the Central Universities (NO. 2017ZY27, 2016ZCQ08, and 2015ZCQ-GX-01).