Abstract
The edge-detection method based on polynomial annihilation that was recently proposed has been applied to locate small damage in structures and demonstrated its effectiveness on beam-like structures. However, significant computational effort involved in this method lengthens the damage detection process, which forbids real-time damage detection. To alleviate this difficulty, in this article, we improve the method suggested by Surace and colleagues by first using the divided difference approach on the identified mode shapes to identify the regions in which jump discontinuities are potentially located and then only applying the polynomial annihilation derivative detector to data points in the identified regions. In this way, the computational burden of this approach is significantly relieved, while the accuracy of damage location is still maintained. The improved method has been validated by numerical simulations on a complex cable-stayed bridge model. This approach does not require baseline response data of structures.
Introduction
Cable-stayed bridges have become the first priority to cover long spans of gaps, such as a river, valley, canal or another road. They are attractive due to the capability of spanning long distances, limiting and eliminating piers (self-anchoring), slender decks, flexibility in the construction schedule, and minimal environmental impact. In addition, they fascinate architects by providing a chance to apply their rich imaginations to create special esthetic shapes. As structures for transporting hundreds or even thousands of people every day, it is important to maintain the integrity and safety of cable-stayed bridges. During their service lives, deterioration/damage may accumulate in bridges due to material aging, environmental corrosion, and extreme loadings, such as strong winds and earthquakes. Failure to detect small damage occurring in the bridge at an early stage may eventually endanger the integrity of the overall bridge.
Along with the increase in the complexity and size of structures, global vibration-based structural damage detection approaches exhibit considerable promise (Doebling et al., 1998; Salawu, 1997; Sohn et al., 2003). Their major findings are presented as follows. Various modal parameters, such as natural frequencies and mode shapes as well as the derivatives of mode shapes, have been widely applied for damage detection. In particular, changes in natural frequencies have limitations on locating damage, because structural damage in different locations may cause the same change in natural frequencies. Therefore, mode shapes were considered as better indicators of damage due to their ability of containing spatial information on damage sites. Later on, curvatures of mode shapes have been found to be more sensitive to damage, especially to minor damage (Roy and Ray-Chaudhuri, 2013; Whalen, 2008). By computing the difference in the curvature of a mode shape before and after damage, peaks associated with damage sites stood out (Pandey et al., 1991). However, this approach requires that mode shapes of the intact structure are available, which is unfortunately not the case for most existing structures. Then some effort has been made to try to overcome this problem. Ratcliffe (2000) proposed a “gapped-smoothing technique,” which was based on a damage index that is a function of curvatures of mode shapes in the current (damaged) state. This approach was susceptible to measurement noises in a mode shape (Kim et al., 2006). Gokdag and Kopmaz (2009) proposed an approach to obtain the information of the undamaged state from damaged mode shapes, and thus, undamaged mode shapes were not required. In particular, recent applications of edge-detection methods in structural damage detection have demonstrated that the edge-detection methods may offer a potential to detect small damage at an early stage (Surace et al., 2013).
Edge detection has been widely used in the fields of image processing and computer vision. In those applications, the purpose of edge detection is to localize variations of the gray level of an image and to identify the corresponding physical phenomena behind the variations (Djemel and Salvatore, 1998). Most edge detectors use the local pixel intensity gradient, which can be obtained by calculating the difference of the convolution of weighted matrix (called local gradient mask) and the image. These detectors include well-known edge detectors, such as Sobel, Roberts, Prewitt, Robinson, Kirsch, Frei-Chen, and Marr-Hildreth (Dusan and Damjan, 2007). However, most of these edge detectors are based on the first-order accuracy of a function. And the performance evaluation is usually judged subjectively according to different users and different requirements for the same image (Kang and Wang, 2007). Even the most popular Canny edge detector (Canny, 1986) still involves the use of several parameters that have to be adjusted in a specific implementation (Saxena, 2008).
Recently, the polynomial annihilation edge-detection method was proposed (Archibald et al., 2005; Saxena, 2008). It is to detect jump discontinuities in a piecewise smooth function and in its derivatives. Different from traditional edge detectors, the polynomial annihilation edge-detection method performs a high-order reconstruction of the jump function, which is an approximation to the difference in function values between consecutive data points. The high-order design of the polynomial annihilation edge-detection method makes it more robust when the underlying structure of the piecewise smooth function has some variability. Moreover, this method does not require the data points to distribute uniformly. In the study by Surace et al. (2013), this method has been applied in locating cracks in beam-like structures when only a few post-damage mode shapes are available. Its fundamental idea is to detect damage by locating jump discontinuities in the first derivative of mode shapes. By combining polynomial annihilation derivative detector, minmod limiter (Gelb and Tadmor, 2006), and stencil shifting technique (Archibald et al., 2008), cracks can be localized through removing oscillations near the discontinuities and clearly locating jump discontinuities. Although cracks can be accurately located, due to the use of the stencil shifting technique, the polynomial annihilation derivative detector has to be run repeatedly many times for each reconstruction data point, which lengthens the damage detection process. In addition, the approach was only applied to simple structural components, such as beams. In this study, the approach will be improved by first identifying the regions in which jump discontinuities are located and then employing the polynomial annihilation derivative detector on only the data points in the identified regions. In this way, the computational effort will be significantly reduced and the damage detection process will be shortened, while the detection accuracy of damage is maintained. Then, this improved approach will be applied to locate damage on a complex structure, for example, a cable-stayed bridge model. The improved approach will make a big difference in the data-processing time for damage detection of complex structures where a lot of data points may be involved.
Review of polynomial annihilation derivative detector
The polynomial annihilation derivative edge detector is briefly reviewed in this section. It is developed based on the polynomial annihilation edge-detection method. The method is to locate jumps in the derivatives of a function (Archibald et al., 2005; Saxena, 2008).
Consider a piecewise smooth function f ∈ Cγ−1, γ∈N = {1, 2,…}, known only on the set of discrete data points, S = {x1, x2,…, xN} ⊂ [a, b]. Assume that f(x) and all its derivatives up to f(γ−1)(x) are continuous in [a, b] and the jump discontinuity first appears in f(γ)(x). For any data point x in the domain, f(γ)(x) has well-defined one-sided limits. Denote the set of jump discontinuities in the γth derivative as Jγ
where f(γ)(x+) and f(γ)(x−) are the right and left side limits of f(γ)(x). The local jump function for the γth derivative is defined as
Therefore, for x ∉ Jγ, [f(γ)](x) = 0, while for x = ξ ∈ Jγ, [f(γ)](x) = [f(γ)](ξ). Consider any reconstruction data point x ∈ (a, b), typically chosen as the mid-points of two adjacent data points, and let m > γ be a positive integer, then the selected local stencil can be expressed as
where Sx is the set of sampled data points around x, so that Sx is composed of m + γ + 2 data points around x. Define
Note that the stencil Sx should be carefully chosen so that
The coefficients cj(x) can be determined by solving the following linear system
Subject to the constraints
where h(x) is defined as
Here, pl, l = 0,…, m is the basis of ∏ m , the space of polynomials up to degree m. The scaling factor qm,γ(x) in equation (5) can be defined as
Like other high-order edge-detection methods, the polynomial annihilation derivative detector has the following problems: when choosing a small m > γ, steep gradients may be falsely detected as jump discontinuities and when choosing a large m, oscillations in the vicinity of a jump discontinuity are inevitable. Therefore, the minmod limiter was introduced in Gelb and Tadmor (2006). It is defined as
where {f1, f2,…, fm} is a given finite set of functions defined in the interval [a, b]. However, the minmod limiter sometimes cannot remove all the oscillations, because at some points near the discontinuity, all oscillations may be in the same direction. To solve this problem, a stencil shifting technique was introduced in Archibald et al. (2008). The stencil shifting technique shifts the original stencil Sx one point left and one point right, respectively. Therefore, three stencils will be obtained, which can be denoted as
When applying the above derivative detector to detect a crack in a beam, it is assumed that there exists a jump discontinuity in the first derivative of the mode shape caused by a crack. Therefore, the first derivative of a mode shape of the beam with a crack can be considered as a piecewise smooth function. To detect crack, it is to identify the discontinuity in the first derivative of a mode shape. To illustrate this, an example in Surace et al. (2013) is considered here. The example will also be used in the next section to compare the efficiency of the improved method and previous methods. For detailed information on the example, the readers are referred to Surace et al. (2013).
Example 1. Consider a cracked cantilevered beam with the following mechanical and geometrical characteristics: Young’s Modulus E = 2.06 × 1011 N/m2, density ρ = 7850 kg/m3, length L = 1 m and squared cross section with b = 0.04 m, as shown in Figure 1. It is assumed that the crack maintains open. For both the intact and damage cases, modal analysis is performed to obtain the first few mode shapes. The first mode shape and its derivative are presented in Figure 2. Although there exists a jump discontinuity in Figure 2(b), it is very difficult to locate this discontinuity by visual inspection. This is why the derivative detector is introduced here. The derivative detector is applied to the first derivative of the first mode shape here, and the obtained crack detection results are presented in Figure 3.

Finite element model of the cracked beam used in Example 1.

First mode shape and its first derivative of the undamaged and damaged beams used in Example 1: (a) first mode shape and (b) first derivative of the first mode shape.

Damage detection results of Example 1 using (a) polynomial annihilation derivative detector with a fixed stencil
Figure 3(a) presents the results when a minmod limiter is used, while Figure 3(b) presents the results when both a minmod limiter and the stencil shifting technique are used. By comparing Figure 3(a) and (b), we can observe that combining the polynomial annihilation derivative detector with a minmod limiter and the stencil shifting technique can eliminate the oscillation in the vicinity of the damage, facilitating in locating damage. However, from our previous analysis, the drawback of using both a minmod limiter and the stencil shifting technique in the polynomial annihilation derivative detector is that for every reconstruction data point, the derivative detector has to be applied repeatedly for each value of m and each shifted stencil. That is, significant computational effort is demanded, especially when a complex structure with a great number of data points is dealt with.
For comparison, the mode shape curvature-based approach is also used to locate the crack, and the results are displayed as the dashed green graphs in Figure 3. Damage can be indicated by the peaks in the difference in the mode shape curvature before and after damage (Pandey et al., 1991). Unfortunately, in this case, a peak occurs where no crack takes place. In addition, this approach requires the mode shapes before damage. By contrast, the derivative detector-based approach only requires the mode shapes in the current state.
Improvement of polynomial annihilation derivative edge detector
Since a crack in a structure is implicated as a jump discontinuity in the first derivative of a mode shape (Surace et al., 2013), we will only focus on the case where γ = 1 in the following. To satisfy m > γ, the value of m should be equal to or greater than 2, and the size of a stencil should satisfy the condition of m + γ + 2 ≥ 5. Previous research (Surace et al., 2013) has suggested that the polynomial annihilation derivative detector combined with the minmod limiter and stencil shifting technique successfully localizes all jump discontinuities. However, as mentioned above, this procedure has to be applied repeatedly for every reconstruction data point, which demands significant computational effort. To alleviate this difficulty, the authors propose to first identify the region in which a jump is located, and then for points outside the region, Lm, γ f(x) are simply set to 0, because an ideal Lm, γ f(x) should have the following property
where Ix is the smallest closed interval such that Sx ⊂ Ix, with Sx defined in equation (3).
In Hu and Shu (1999), the divided difference of a function f(x): [a, b] → R is used to measure its smoothness. It is well known that the first-order divided difference is defined as f[xi] = f(xi), ∀i = 1,…, Ns, where Ns is the size of a stencil. The kth-order divided difference of f(x) is recursively defined in terms of the (k−1)th-order divided differences as
According to Hu and Shu (1999), the divided difference of f(x) satisfies the following theorem.
Theorem 3.1. If f(x) is smooth inside a stencil Sx = {x1,…, xm+γ+2}, then
where ξ ∈ (x1,…, xm+γ+2). If f(x) is discontinuous at some point inside the stencil, then
where the definition of h is given in equation (8).
According to Theorem 3.1, if f(x) is smooth inside a stencil Sx, then f(x) has derivatives of all orders inside Sx. Thus, the first derivative of f(x) is also smooth inside Sx. In this case, equation (13) is valid, and f[Sx] will approach to 0 with the increase in m. Otherwise, if Sx encloses a discontinuity in the first derivative, then, by Theorem 3.1, equation (14) becomes valid. Since h is always greater than 0, we have f[Sx] = O(1/hm+γ+2) > 0. Based on the above criteria, the divided difference can be used to identify whether a stencil encloses a discontinuity in the first derivative of f(x) or not. The smaller the divided difference is, the smoother the first derivative will be. And this will become more evident with the increase in m. This can be illustrated by applying the divided difference operation to the function in Example 1. The divided difference results are presented in Figure 4(a) and (b). This figure shows that there is a dramatic difference in the divided difference between the region in which a jump discontinuity is located and the smooth regions. Based on this observation, we proposed the following algorithm (Algorithm 3.1) to identify the region in which a jump discontinuity is located.

Apply Algorithms 3.1 and 3.2 to Example 1 to identify the potential discontinuous regions of the first derivative of f(x). The red diamonds indicate the potential jump discontinuities. The two thresholds T1 and T2 in (a) and (b), respectively, are obtained through Algorithm 3.2: (a) |f[Sy]|, #Sy = m + γ + 2, m = 2, (b) |f[Sy]|, #Sy = m + γ + 2, m = 5 and (c) the list L obtained by using Algorithm 3.1.
Algorithm 3.1. Given a sampled data set
1: Create an empty matrix M ∈ ℛ2×Nr and an empty list L ∈ ℛ1×Nr, set two thresholds T1 and T2;
2: For i = 1, 2:
2.1: For j = 1,…, Nr:
2.1.1: If i == 1, then set m = 2, else set m = 5;
2.1.2: Set y = yj, choose three stencils,
3: For j = 1,…, Nr:
3.1: If M(1, j) > T1, then set M(1, j) = 1, else set M(1, j) = 0.
3.2: If M(2, j) > T2, then set M(2, j) = 1, else set M(1, j) = 0.
4: For j = 1,…, Nr:
4.1: If M(1, j) & M(2, j) == 1, then set L(j) = 1, else set L(j) = 0.
In this procedure, only two representative cases (m = 2 and m = 5) are considered. When using Algorithm 3.1, it is critical to properly set the two thresholds. The approach proposed in Gonzalez et al. (2009) is adopted in Algorithm 3.2 with some modifications. Because this approach should be executed following Step 2 in Algorithm 3.1, the matrix M is assumed to have been created at this stage.
Algorithm 3.2. This algorithm produces two thresholds T1 and T2, which can be used in Step 3 in Algorithm 3.1. This algorithm is detailed as follows:
1: For i = 1, 2:
1.1: Set T be the median value of the ith row of M;
1.2: Separate the ith row of M into two groups: group G1 contains all the elements with values greater or equal than T, while group G2 contains all the elements with values less than T;
1.3: Compute the average value of µ1 and μ2 for groups G1 and G2;
1.4: Set
1.5: Repeat Step 1.2 to Step 1.4 until the difference in T in consecutive iterations is smaller than a tolerance τ0;
1.6: If i = 1, then set T1 = α × T, else set T2 = α × T, where α is a predefined parameter in [0, 1].
Using Algorithms 3.1 and 3.2, data points in the vicinity of a jump discontinuity can be identified. Figure 4(c) plots the list L obtained by applying the two algorithms to Example 1. Based on the obtained list L, the following detection procedure will be followed: for any reconstruction data point yj ∈ Y, if its corresponding Lj = 0, then Lm, γ ,Syf(yj) is set to 0; for the regions with Lj = 1, the polynomial annihilation derivative detector with multiple stencils (PAMS) will be employed to detect jump discontinuities. Herein, the only parameter to determine is α in Algorithm 3.2, when m is assumed to be 2, 3, 4, and 5, and τ0 is fixed as 10−8. Figure 5(a) presents the result by applying the above procedure to Example 1 with α = 0.6. By comparing Figures 5 and 3, the improved method can effectively remove oscillations near jump discontinuities and pinpoint the discontinuity in the function and then the crack location. The accuracy of the improved method is the same as the previous method, PAMS, where the PAMS has to be applied to each construction data point in Y repeatedly.

Damage detection results using the improved method and comparison of time consumption among methods: (a) damage detection result based on the list L, red diamonds indicate the potential discontinuities (
To demonstrate the efficiency of the improved approach, Figure 5(b) presents the time consumption of each method with the increase in the size of reconstruction data set (from 25 to 50). It can be seen that the improved method (DivDiffPAMS) consumes much less amount of time than PAMS and consumes similar amount of time to the polynomial annihilation derivative detector with a fixed stencil (PAFS), when the number of data points is 25. When the number of data points is 50, the advantage of the improved method in shortening the data-processing time becomes more obvious. In this case, the improved method even consumes less amount of time than PAFS. In summary, the improved method significantly reduces the computational effort and time, while maintaining the same damage detection accuracy as PAMS and possessing better accuracy than PAFS.
Numerical simulation
In this section, the improved approach is applied to a two-span continuous cable-stayed bridge model to validate its effectiveness. A finite element model of this bridge model is developed using ANSYS, as shown in Figure 6. The left and right spans are 6.5 and 2 m long, respectively. This bridge model is 0.8 m wide. The two piers are 0.4 m high, and the two towers are 2.6 m high from the deck. It has two panels of harped cables with 17 cables on each panel. The bridge deck system on each span is a steel sheet with a thickness of 2 mm supported by a steel grid. The steel grid is composed of two longitudinal girders with an H-shaped cross section and transversal beams with a rectangular cross section at an interval of 0.25 m in the longitudinal direction. Each girder is rigidly connected to a pier and is supported by rollers at the two ends. The piers and towers are assumed to be constructed by steel with an H-shaped cross section. Cables are assumed to be constructed by steel strands. Modal analysis is first performed to obtain the first four natural frequencies and mode shapes of the bridge model, as shown in Figure 7.

A finite element model of the cable-stayed bridge model.

First four mode shapes and natural frequencies (a) first mode shape (first bending mode, 9.64 Hz), (b) second mode shape (first torsion mode, 13.07 Hz), (c) third mode shape (second bending mode, 24.57 Hz), and (d) fourth mode shape (second torsion mode, 29.31 Hz).
Damage is simulated by adding steel plates underneath the bottom flange of the girder on the front panel. Two damage scenarios are simulated here. In Damage Case 1, a steel plate is added between Nodes 28 and 29, as shown in Figure 8(a). The width of the added plate is the same as the width of the girder, 0.1 m; its length is equal to the distance between two beams, 0.25 m, and its thickness is the same as the thickness of the flange of the girder, 0.008 m. In Damage Case 2, one more steel plate with the same dimensions is added between Nodes 17 and 18, as shown in Figure 8(b), while keeping the added steel plated in Damage Case 1. Please note that the thickness of the added plate is exaggerated in the drawing of Figure 8. Pre-analysis has been conducted to ensure that the simulated damage is away from the node points of the first few mode shapes.

Two simulated damage scenarios: (a) Damage Case 1 and (b) Damage Case 2.
Damage detection results
For each damage case, the improved method is applied on mode shapes acquired on nodes from 11 to 38 by setting m = 2, 3, 4, and 5, and α = 0.6. Reconstruction points are uniformly distributed, and the interval between two reconstruction points, that is, the value of h(x) in equation (8), is set to 0.5, so there are totally 55 reconstruction points. For comparison, the PAFS and PAMS are also employed. As in the previous section, we use “PAFS” to represent polynomial annihilation derivative detector with a fixed stencil, “PAMS” to represent polynomial annihilation derivative detector with multiple stencils, and “DivDiffPAMS” to represent the improved method developed in this study.
For Damage Case 1, because similar results are obtained when using a different mode shape, only representative results using the first and second mode shapes are presented here. Figure 9 presents the first derivatives of the first and second mode shapes in the intact and damaged states. It can be seen that damage does not appreciably change the mode shape derivatives, and thus, damage cannot be identified from them visually. To apply the improved method, Algorithms 3.1 and 3.2 are first applied to all mode shape components to find the potential region in which jump discontinuity is located. The identified region is indicated by red diamonds in Figures 10(c) and 11(c). Then, the polynomial annihilation derivative detector is only applied to the identified region to locate the damage. In Figures 10(c) and 11(c), a spike occurs to node 28, which indicates that damage occurs between nodes 28 and 29 (specified in Damage Case 1). Please note that the horizontal axis of Figures 10 and 11 represents data at the interval of 0.5, since the interval between two adjacent reconstruction points is 0.5. For comparison, the results obtained by the mode shape curvature-based method (dashed green lines) are not as good as the improved method.

First derivatives of mode shapes before and after damage: (a) first derivative of first mode shape and (b) first derivative of the second mode shape.

Damage detection results for Damage Case 1 using (a) PAFS, (b) PAMS, and (c) DivDiffPAMS when the first mode shape is used.

Damage detection results for Damage Case 1 using (a) PAFS, (b) PAMS, and (c) DivDiffPAMS when the second mode shape is used.
The damage detection results using the PAFS and PAMS when the first mode shape is used are presented in Figure 10(a) and (b); the results when the second mode shape is used are presented in Figure 11(a) and (b). By comparing Figure 10(a) with (c), and Figure 11(a) with (c), it can be seen that all the oscillations near the jump discontinuity when using the PAFS have been removed by the improved approach. Although equally good results are obtained when using PAMS, as shown in Figures 10(b) and 11(b), it takes much more time to run the PAMS than the improved method.
For Damage Case 2, representative damage detection results using the third mode shape are presented in Figure 12. The two damage sites are indicated by two spikes. In Figure 12, spikes occur at nodes 18 and 28, which indicates multiple damage sites (between nodes 17 and 18, and between nodes 28 and 29). It demonstrates that the improved method is also suitable for identifying multiple damage sites. Although both the PAMS and the improved method deliver good results, the polynomial annihilation derivative detector is only run on a less number of data points when applying the improved method and thus shortens the data-processing time.

Damage detection results for Damage Case 2 using (a) PAMS and (b) DivDiffPAMS when the third mode shape is used.
Conclusion
The polynomial annihilation derivative detector has been proven to be an efficient method to locate damage in beam-like structures when only a few post-damage mode shapes are available. However, it has to be applied repeatedly a number of times to each reconstruction data point when using a stencil shifting technique to remove oscillations in the vicinity of a crack in the damage indicator, which is time-demanding. In this study, the polynomial annihilation edge-detection method is improved to speed up the detection process. It is achieved by first identifying the regions in which jump discontinuities are located and then only applying the polynomial annihilation edge detection to the data points in the identified regions. In this way, the computational burden can be significantly relieved, and the duration for the detecting process will be shortened, while the detection accuracy of damage is maintained. Numerical simulations on a complex cable-stayed bridge model have demonstrated the efficiency of the improved method. This method only requires a few post-damage mode shapes related to lower natural frequencies.
Since the improved approach is based on the change in mode shapes caused by damage, if damage happens to occur at the node points of a specific mode shape, there may be difficulty in detecting this type of damage. Fortunately, the node points of a mode shape may not be the node points of another mode shape. As long as the proposed method is applied to more than one mode shape, damage may still be detected. In this case, the damage detection results from different mode shapes may vary, and the difference may lie in the damage sites at the node points of a mode shape. In addition, this method may be sensitive to the density of measurement points, which will be investigated in future research.
Footnotes
Acknowledgements
The first author (G.Y.) would like to appreciate USDOT through University Transportation Center for their financial support.
Declaration of conflicting interests
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The first author (G.Y.) also would like to thank University of Texas at El Paso (her previous institute) for their consistent support, as most of the research work was conducted there.
