Noise Reduction for Improvement of Ultrasonic Monitoring Using Coda Wave Interferometry on a Real Bridge

Reinforced concrete bridges are iconic parts of modern infrastructure. They are designed for a minimum service life of 100 years. However, environmental factors and/or inappropriate use might cause overload and accelerate the deterioration of bridges. In extreme cases, bridges could collapse when necessary maintenance lacks. Thus, the permanent monitoring for structure health assessment has been proposed, which is the aim of structural health monitoring (SHM). Studies in laboratories have shown that ultrasonic (US) coda wave interferometry (CWI) using diffuse waves has high sensitivity and reliability to detect subtle changes in concrete structures. The creation of micro-cracks might be recognized at an early stage. Moreover, large-volume structures can be monitored with a relatively small number of US transducers. However, it is still a challenge to implement the CWI method in real SHM practical applications in an outdoor environment because of the complex external factors, such as various noise sources that interfere with the recorded signals. In this paper, monitoring data from a 36-m long bridge girder in Gliwice, Poland, instrumented with embedded US transducers, thermistors, and vibrating wire strain gauges, is presented. Noise estimation and reduction methods are discussed, and the influence of traffic, as well as temperature variation, are studied. As a result, the relative velocity variation of US waves following the temperature change with a very high precision of 10% is shown, and a good bridge health condition is inferred. The influence of lightweight real traffic is negligible. The study verified the feasibility of the implementation of the CWI method on real bridge structures.


Introduction
For many decades, concrete has been widely used for civil infrastructure construction (i.e., bridges, highways, foundations, etc.) because it is strong, durable, and resistive. The bridge construction market is growing rapidly to meet the demand of economic development. Due to the increasing freight traffic, environmental factors and/or irregular use, the service life time of bridges can be reduced significantly. Accordingly, the timely necessary maintenance of bridges to guarantee their service life is receiving more attention. Damage detection at a very early stage will undoubtedly reduce the maintenance costs and prevent the failure or complete collapse of bridges. Both temporary non-destructive evaluation (NDE) and permanent long-term structural health monitoring (SHM) play a more critical role. Researchers have been working on structural health monitoring since the 1970s; however, there is still a gap between research and industrial applications [9].
A wide variety of sensors have been developed for bridge health monitoring. Fiber optic sensors can provide distributed strain information [7]a sw e l la st h ei n f o r m a t i o n about temperature variation [4,31]. Local strain can also be measured by classic strain gauges (SG) [5,6] or advanced high-sensitivity wireless strain sensors [21]. Accelerometers, which measure the vibration acceleration data, are another type of common sensor used for bridge SHM applications [19,20,22]. In recent research, along with the development of technologies, smartphones can also be involved in SHM. Some sensors integrated inside the smartphone may contribute valuable information for SHM [11,30]. Moreover, most modern smartphones' hardware capacity and perfor-mance are sufficient for data analysis, if the corresponding software is adapted [17]. However, most of these sensors are installed or glued at the surface of the structures. The influence of the external environment can not be ignored. The application of SG is a point-based measurement; it provides a local strain variation measurement around a fixed position; defects cannot be detected when they appear in other areas away from the SG. If a long girder should be monitored, lots of SGs are needed. Optical fibers usually are fragile, and the data acquisition is complicated and expensive. To overcome these problems, new techniques based on using sonic methods are proposed.
The sonic technique is a promising NDE technology as the wave velocity can be directly related to the material's physical and mechanical properties. The Ultrasonic Pulse Velocity (UPV) test is an in-situ test to evaluate the quality of concrete by measuring the wave propagation velocity [10,27,36] in a concrete body. Here, a higher velocity indicates better quality [23]. Though this method has a limitation in detecting weak velocity variations, it is still one of the most successful techniques. Many commercial products are based on this technique and are widely used for industrial applications. However, this technique can usually not be applied for long-term SHM. Acoustic Emission (AE) is a passive sonic method to identify defects when they develop. It can detect damage at a very early stage [16]. This technique has been successfully applied for real bridge SHM [42,45]. However, as the AE signals are relatively weak, they are often strongly disturbed by the noise in an outdoor environment, such as the traffic-induced noise. Moreover, AE requires a high data acquisition rate for data collection; therefore, a large amount of data is recorded during AE monitoring, and costly data analysis is required.
The research about SHM on bridges focuses mostly on sensor development and development of algorithms for data analysis [14]. In this paper, a new type of embedded US sensors is presented, and a recent method, named the Coda Wave Interferometry (CWI), is introduced to monitor a real bridge in Poland. Embedded US transducers were installed inside the structure in order to reduce the influence coming from the external bridge environment. This way, the internal changes in the structure will then be better monitored. The CWI method was initially proposed for seismology research [1,2]. Since the year 2000, it was progressively introduced into the NDE domain [39]. This method extracts useful information by analysing the later parts of diffuse waves in the end part of recorded signals. Existing studies have shown its extremely high sensitivity in detecting weak changes in the wave propagation medium, which affects the coda signal through a perturbation of its waveform and/or a shift in the arrival time of certain wave train. However, it is still a challenge to implement the CWI method on real concrete objects such as bridges for long-term monitoring. Other than in the laboratory, the outdoor environment poses multiple influences that make the analysis more difficult. This paper describes how the CWI method is applied at a bridge, which is equipped with a new type of embedded US transducer. Noise estimation and reduction are investigated, and velocity variations, induced by temperature changes and strain, as well as the influence of the traffic passing the bridge, are discussed.

Thermally-and Strain-Induced Velocity Variation
The condition of civil infrastructures is mainly influenced by two main factors: temperature variation and load effect. Wave velocity is strongly related to the elastic properties of the propagation medium. The linear relationship between wave propagation velocity and elastic strain when a limited load is acting on the non-linear elastic material (mainly concrete) object has been confirmed in [3,43,48]. When the strain variation exceeds a certain level, and damage such as microcracking appeared, the linear relationship is broken. A temperature variation may lead to a geometry and/or elastic properties change of the structure being monitored. Plenty of studies have shown the influence of temperature on the elastic modulus and strength of concrete material [13,26,29]. In other words, the temperature variation can indirectly create internal strain changes, leading to a velocity variation. Depending on different situations, the temperature effect could be an unwanted factor or a valuable source. In most cases, the temperature effect needs to be eliminated, and a thermal bias control technique is presented in [49]. On the other hand, temperature-induced strain data could also be used for damage assessment [46]. Grêt et al. used thermally-induced velocity variation to detect microcracking in [15]. The damage level of a concrete specimen was also identified by Sun and Zhu in [41] using the relation between wave propagation velocity and temperature. In this paper, temperature variation is the most crucial variable, and the wave propagation velocity variation is mainly thermally induced.

Coda Wave Interferometry
CWI method is a promising technique to monitor weak wave velocity variations in a heterogeneous medium. When US measurement is carried out in a frequency range exceeding 50 kHz, multiple scattering is caused by the heterogeneity  [35]. Waves travel along much longer paths through a wider region compared to the direct wave. The impacts of the weak changes in the medium are accumulated, which ensures the high sensitivity of the CWI method in detecting subtle changes. CWI has been successfully applied in laboratory environments to detect stress variation [40,48], temperature fluctuations [34,39] and cracks at a very early stage [3,43] in small-scale specimens. Furthermore, detection of stress-induced velocity change at a large-scale specimen under field conditions [44] and at a real bridge [40] was also achieved. Figure 1 shows two US signals recorded before and after a weak change in the propagation medium. As one can see, the first arrivals of the two signals are identical; however, a slight perturbation can be found in the later arrivals. The fundamental of the CWI method is to extract two useful features, relative velocity change (dV/V) and correlation coefficient (CC), from two US signals measured in different states. The CC measures the similarity between two signals. It can be used for imaging purpose [24,47,50], or to monitor stress change [18,33]. Various CWI processing approaches have been developed, such as the Doublet techniques [39], Taylor series expansion method [28], and the stretching method [32]. In this paper, the stretching method is chosen for CWI analysis as it provides the most accurate results among these three approaches [49]. The time shift is considered as a time dilation or compression by a factor α, indicating a velocity decrease or increase. u u (t) and u p (t) represent the signals recorded in two different states. The α, which maximizes the cross-correlation (Eq. 1) of the two signals within the time window [− T, T], is chosen as the relative velocity change. Window selection for this experiment is presented in Sect. 4.1.

Gliwice Bridge Gm1-2W
The Gm1-2W bridge is a post-tensioned reinforced concrete bridge constructed in 2014 as an extension of Highway 902, connecting Katowice city and Gliwice city in Poland. The bridge's total length is 552 m, and it consists of ten interior 48 m long spans and two external 36 m long spans. It is supported by piers, piles, and two girders 8 m apart, as shown in Fig. 2. The monitored part of this bridge was one of the 36 m long, 15.1 m wide external spans located at the south-east end of the Gm1-2W bridge. As the middle part of the span undergoes the most substantial stress changes, the centre part of the girder on the west side was selected for the investigation (Fig. 3). Ultrasonic transducers, thermistors, and vibration wire strain gauges were integrated into this area during its construction before the casting in 2014. The experiment presented in this paper was the first formal analysis on this bridge girder using the combination of all these sensors. The bridge girder was monitored continuously for 5 days to study the changes or damage due to traffic and environmental influences. It was a first step to validate the methodology for long term monitoring.

Transducer SO807
"SO807" is a special type of embedded US transducer designed by Acoustic Control Systems, Ltd. exclusively for BAM research [33]. The main body of SO807 is a piezo-ceramic cylinder. It generates US signals in the directions which are perpendicular to its vertical axis. Because of the directivity limit, two transducers installed in parallel provide the ideal measurement setup. Compared to the common piezoelectric bender elements, SO807 can generate a more powerful and stable US signal, which meets the needs of long-term measurement on real structures. The detailed dimensions of SO807 are presented in Fig. 4.A l l the electric connections and cables are inside the hollow interior of SO807, which ensures better protection during installation. This transducer can act both as transmitter and receiver. The central frequency of SO807 is around 62 kHz. The wavelength in concrete related to this frequency usually is two times larger than the size of most aggregates [12]. Thus, waves interact with the heterogeneities. A combination of L-shaped rebar and a PVC tube was used to install the US transducer (Fig. 5). One end of the L-shaped rebar was welded to the reinforcement to stabilize the position of SO807. The depth position of SO807 was given through a PVC tube connected to its end metallic end. This installation method ensures that the transducer will stay at its specific position during the concrete casting and vibration process. Seven US transducers, three thermistors, one rebar strain sensor, and one vibration wire strain gauge were mounted in the monitoring area during the construction before casting. Positions of all the sensors are shown in Fig. 6.Asthe sensors were installed inside the structure, protected from the external influences, they were still functioning well four years after the construction. All seven US transducers were installed in a perpendicular direction to the bridge. The distance for different transmitter-  Table 1. The arrangement allows both long-range measurements parallel to the girder's main axis and short-range measurements. The most extended range presented in this paper is around 2.19 m, and the shortest range is 0.77 m.

Data Acquisition System
The data acquisition system (Fig. 7) consists of a laptop with a custom-built control software to store data and configure the measurement parameters such as measurement interval and transmitter-receiver combination order; a digital-analogue data acquisition system to convert the signal and to trigger the excitation; an amplifier to amplify the excitation signal; a pre-amplifier with analogue filter to pre-process the signal from the receiving transducer; a multiplexer to switch between different transducer combinations and a power supply. To avoid theft and for safety reasons, all devices were locked in a protection box 4 m above the ground (Fig. 8). This data acquisition system allows a continuous measurement with a sampling frequency of 1 MHz. In this experiment, the time required for one round of acquisition (twelve transducer combinations) was 33 s. The bridge was continuously monitored for 5 days, from October 26th to 30th 2018. More than 12,000 signals were recorded for each transducer combination. Temperature and strain data were collected by a Campbell data logger every three minutes.

Results
The active measurements were performed for 5 days to detect changes in the structure. It was assuming that the humidity and corrosion conditions remained the same. The major change in the structure was then the temperature. The system acquired raw ultrasonic signals, including the traffic-induced   It is necessary to analyse the raw signals and convert them into meaningful information indicating the changes in the structure. As a first thing, it is essential to distinguish the noise from the acquired signals and then reduce it to a minimum level.  Figure 9 shows a raw US signal recorded by transducer combination S06E08 at midnight on October 27th 2018 as an example. One hundred samples (0.1 ms) were recorded before the signal transmission ('pre-trigger') to determine the noise level and offset. Due to electromagnetic interference, crosstalk was recorded at 0.1 ms when the excitation signal was sent. A pre-processing procedure is required to remove offset and crosstalk and to reduce the noise. In a laboratory environment, noise can be filtered in most cases by applying a simple band-pass filter (e.g., passband [10 kHz, 150 kHz]). Details of this procedure are presented in [32]. The evaluation can be limited to the part of the signal containing significant US information. Thus, a suitable window should be selected to improve the efficiency of data processing. In this study, window [600, 1600] ([0.5 ms, 1.5 ms]) is chosen for data analysis and signal-to-noise ratio (SNR) level estimation. Assuming that the same noise (the first 100 samples) runs through the entire signal, SNR could be calculated

Noise Estimation
The SNR level of all signals recorded by transducer combination S01E02 for 5 days is shown in Fig. 10 as an example. A significant difference is observed. During the day (from 8 a.m. to 6 p.m.), the average SNR level was around 51.18 dB while it decreased to 42.84 at night. The noise level increased a lot at night. This phenomenon also appears in the measurements by the other transducer combinations. The authors assume that the power supply for the street lights on the bridge caused this difference and created additional electromagnetic noise as soon as the lights were turned on. Hence, the main noise in the measurement on this bridge came from the high voltage electric power supply passing through the lamp cables. This effect had never been observed in the laboratory environment. The average SNR level of all transducer combinations is presented in Fig. 11. This result proves that the noise level became higher during the night. Unfortunately, the SNR levels of transducer combinations S04E07, S05E08, and S05E07 are unexpectedly low. For example, the SNR level of S05E07 is 8 dB, which means the signal's power was only 6.3 times stronger than that of the noise. The weak signals affects the CWI evaluation. By examining the signals recorded by all the transducer combinations, the noises are almost identical over the whole monitoring area; however, the signal energy was enormously depending on the deviation angle between two US transducers and the coupling condition between the transducers and the concrete. The most intuitive way to study the noise is to analyze it in both time and frequency domain directly. As shown in Fig. 9, the last 1 ms (1000 samples) of the raw data is basically noise. Sixteen signals (eight taken during the day and eight during the night) recorded by transducer combination  Fig. 12. The amplitudes of the noises recorded during the day fluctuated in the range ±20, while they went up to ±150 at night. All the frequency components of the noises recorded during the day were kept at a low level (less than 1.75). The main frequencies were distributed around 62 kHz, which is the primary frequency of the SO807 transducer. According to the amplitude spectrum of the noises at night, the noises contained two main frequency components: 40 kHz (FFT amplitude can exceed 10) and 67 kHz (FFT amplitude can exceed 40). Another way to estimate the noise components is to compare the noisy signal recorded at night with a reference signal recorded during the day. The amplitude spectrum differences between forty signals recorded at night and one reference signal recorded during the day by transducer combination S05E07 are shown in Fig. 13.T h e frequency changed significantly around 40 kHz and 67 kHz, which is in accordance with the previous result.

Noise Reduction and Filter Design
Three different filters are presented in this section. The basic necessary filter applied to the signal is a normal band-pass Fig. 13 The difference of amplitude spectrum (fast Fourier transform) between forty signals measured at night and a reference signal measured during the day by transducer combination S05E07 filter. It can filter the common low-frequency noises, which are the main reason for the signal 'trend' problem (e.g., raw data in Fig. 15). In this study, [30 kHz, 100 kHz] was chosen as the normal frequency band. The Butterworth filter, which provides a flat frequency response, was used to design the band-pass filter. Butterworth filter has a slow roll-off; thus, the order of the filter was chosen to six, indicating a response decrease of − 36 dB per octave. As mentioned in Sect. 3    on this frequency. The order of this filter was also fixed to six. Another idea was to eliminate the noise by its specific frequencies. Thus, a double-notch filter was proposed. This filter was a convolution of the normal band-pass filter and two IIR band-reject filters. The band-reject filter was provided by Python SciPy library. As one of the main frequencies of the noise (67 kHz) was close to the central frequency of the SO807 transducer (62 kHz), the bandwidth of the filter should be as narrow as possible. Thus, sixty was chosen as the quality factor, which characterizes the notch filter − 3 dB bandwidth. This filter passed most of the frequencies but attenuated the components of 40 kHz and 67 kHz to a very low level. Frequency responses in dB of these filters are presented in Fig. 14. Twenty signals recorded by transducer combination S05E07 on the 27th of October at midnight are chosen as examples. All the signals were recorded within ten minutes. Thus, the temperature variation was negligible. Theoretically, all the signals should be almost identical. Signals filtered by these different types of filters are shown in Fig. 15 along with the raw data. The normal band-pass filter successfully removed the signal trend, which was mainly caused by low-frequency noise. However, the differences between these signals are still evident that there exists a tiny irregular amplitude variation and time shift. This will affect the calculation of CC and dV/V. The narrow band-pass filter filtered the noise effectively. All the processed signals are approximately identical. As this filter focuses on the central frequency of SO807 transducer, most of the other frequency components are eliminated. The filtered signals became smoother. However, they contained less information. The similarity of the later arrives was high; in some cases, the stretching method might have difficulties identifying the direction of the waveform shift. The double-notch filter retained most of the original signal's information along with an excellent noise reduction effect.
All these three filters were tested for CWI analysis in Sect. 4.4.

Temperature and Strain Variation
As shown in Fig. 16, the internal temperatures of the monitored area changed less and lag behind compared to the outside temperature. The vibration wire SG determines the strain by measuring the frequency, which can then be converted to strain by a specific gauge factor. The SG was not temperature compensated. The trend of the strain variation was between outside and inside temperature changes. The outside temperature and stain increased significantly after 29th of October.

CWI Analysis
CC and dV/V measured by transducer combination S01E02 are presented in Fig. 17. As the signals recorded by transducer combination S01E02 had the highest SNR level, the normal band-pass filter's effect is good enough. The doublenotch filter reduced the fluctuation of CC and dV/V slightly. The narrow band-pass filter increased the similarity between all the signals (Fig. 17 blue curve in CC); however, it changed the turning points of both CC and dV/V curves unexpectedly. The authors assume that some parts of the filtered signals have some similarities due to information loss caused by the narrow band-pass filter. For instance, the signals filtered   Fig. 15 shows a similar quasi-periodic signal after 0.6 ms. In this case, the stretching method does incorrect estimations of the time dilation factor. CC and dV/V measured by transducer combination S05E07 are presented in Fig. 18. Signals measured by transducer combination S05E07 have the lowest SNR level. The noises have a noticeable influence on CC and dV/V estimation. Narrow band-pass filter and double-notch filter did improve the CC results; however, the impact of the noise on CC values still exists. This is because of the uncertain and random strong noises which interfere with the similarity between signals. All noises have similar frequency components and energy, but their phases are random. Nevertheless, the behaviour of the double-notch filter is acceptable for processing the dV/V calculation. This is one of the advantages of the CWI stretching method that velocity variation could be estimated even under noisy adverse conditions. By comparing the results using different filters, the double-notch filter  was considered as the most suitable filter in this study for the conditions around the bridge. A Savitzky-Golay filter (window 200 and order 3) was applied to the CC and dV/V curves to ensure the smoothing of the results without distorting its tendency. The filter is also provided by the python Scipy library. It is based on the local polynomial least-squares fitting in the time domain. Details of the filter are presented in [38]. CC and dV/V curves of transducer combination S02E06 filtered by Savitzky-Golay filter are illustrated in Fig. 19 as an example. The result of the filter was remarkable to smooth the result.

Thermally Induced Relative Velocity Change
Thermistor T1, T2, and T3 captured the internal temperature at three different positions along the same line with US transducer numbers 5 and 6 ( Fig. 6). As shown in Fig. 16, these three temperatures have a similar global trend with minor different local variations. An average temperature value T average of these three measurements was calculated to present the temperature change in transducer combination S05E06 measured area. CC and dV/V values of all transducer combinations are shown in Fig. 20 along with the average temperature. CC curves measured by the first six transducer combinations (from S01E02 to S06E08) whose SNR levels at night were higher than 21 dB had a better continuity compared to the other six. The dV/V curves measured by all the transducer combinations have similar trends and vary slightly in the range [− 0.09%, 0.049%] with a high resolution of 10 −4 %. They are in coherence with the average temperature.
As one can see in Fig. 20, the relative velocity variation is mainly thermally induced. The temperature-velocity (T-V) curves of transducer combination S05E06 and S02E05 is displayed in Fig. 21. The T-V curves are divided into two parts for analysis. The first part is the interval from the beginning of the measurement to October 29th. The second part is the interval from 29th of October to the end of the measurement. In this interval, during when the strain and external temperature changed significantly. The slopes of all T-V cures are shown in Fig. 22. There is a slight difference between the slopes of the two parts, which is in accordance with the result in [37]. The slope of the first part of the T-V curve recorded by transducer combination S05E06 is about − 0.019, and the slope of the second part is around − 0.020. However, the slopes of two parts of the T-V curve recorded by transducer combination S02E05 have a more significant difference.
The authors assume two main reasons for the slope difference and curve hysteresis: the T average is considered the internal temperature variation in the area, covered by transducer combination S05E06. Although the internal temperatures in other areas should be similar to T average ,t h e y might still be different. The transducer combination S05E06 measured mainly in the direction perpendicular to the girder's main axis while S02E05 measured in the parallel direction. In general, the strain variation in the direction parallel to the main axis is higher compared to that of the perpendicular direction. When the tensile stress increased faster after the 29th of October, the wave velocity should be reduced more significantly. Consequently, the slope of the second part of T-V curve decreased, as shown in Fig. 22.

Influence of the Real Traffic
Signals measured by transducer combination S01E02 have the highest SNR level, meaning that the noise's influence is the lowest. Thus, the dV/V curve of S01E02 is chosen for illustration (Fig. 23) to study the influence of the traffic passing the bridge. The weak fluctuation is considered as the influence of real traffic. The largest range of dV/V during this one hour is 0.0028%, which is extremely small. The influence of real traffic has relatively little effect on overall structural response [8]. The global dV/V change caused by the traffic is negligible, and it can be reduced or eliminated by the Savitzky-Golay filter.

Discussion and Conclusions
This study has shown the possibility of applying the CWI method at a real bridge under a noisy outdoor environment. Different filters should be designed to reduce the influence of noise in the environment. Three filters are tested and compared in this study. The double-notch filter was the most suitable filter for this bridge environment. Most of the noise impact was successfully minimized to an acceptable level.
The high sensitivity of the CWI method in detecting temperature variation and strain changes inside the bridge has been verified based on the results. The strain is mainly induced by the outside temperature. During 5 days of monitoring, the largest wave propagation velocity change was measured by transducer combination S05E08. It was around − 0.091%, which is still relatively small. Compared with CC, the results of dV/V were more persuasive in a noisy environment. The energy of US signals measured by transducer combination S05E07 at night was only about six times bigger than that of the noise. Despite this, the dV/V measured by this transducer combination still provided a useable result, just like all the other transducer combinations did. All the weak wave velocity changes were estimated with a resolution of 10 −4 %. The dV/V and internal temperature curves are roughly linear. Higher temperature leads to a lower velocity. The internal temperature changes were slower and smaller compared to the external environment temperature. Because of this hysteresis, after the 29th of October, when the strain increased was faster, following the external temperature variation, the internal temperature changes remained in a smaller range. This is the main reason for the slope difference presented in Sect. 4.5.
CWI method detected both internal temperature variation as well as the thermally induced strain change. However, it is still difficult to separate the two influence factors. The external temperature variations strongly influence both internal temperature and strain.
When damage, such as a crack, appears in the transducer monitored area, the loss of wave propagation velocity exceeds 1% [40]. The dV/V measured by all transducer combinations varied in an extremely small range. The slopes of all T-V curves were less than − 0.027%, which is in accordance with existing studies [25,37,41]. According to [41], when the structure's damage level becomes higher, the absolute value of the slope of T-V curve increases. Thus, the authors concluded that the monitored bridge span was in good health condition. The influence of the lightweight traffic on the global evaluation is negligible.
For long-term monitoring, a larger measurement interval, such as one hour, should be chosen to reduce the cost of data storage and calculation. For future implementation and monitoring of other structures, it is suggested to install more thermistors in the monitored area to have a more accurate T-V relationship. Moreover, we can appropriately reduce the number of US transducers in a 4-m-long area, because large structures can be monitored well enough by limited numbers of the embedded US transducers. Especially in combination with other existing SHM techniques, the CWI method can provide substantial improvements for the information gain during long-term SHM at concrete building structures.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.