Chaotic signals have long held promise as a means of excitation in structural health monitoring applications, but methods to process the structural response and infer damage are limited in number and effectiveness. Here, an alternative geometric methodology is presented that is based on measuring the boundary deformation of a system attractor as parameters change. This technique involves sampling the boundaries of two system attractors: one with nominal parameters and one with varied parameters, and then computing boundary transformation vectors (BTVs) between them. These vectors encode information about how the system has changed. This method allows damage level as well as type/location to be simultaneously quantified in simulated structures, and represents a major step toward making chaotic excitation a more practical choice for structural health monitoring.
Being able to monitor the condition of a structure is important in applications ranging from civil infrastructure to rotating machinery. Such structural health monitoring may be concerned with detecting the existence, type, location, and severity of damage [1,2]. Identifying damage in dynamic systems is often accomplished using active, vibration-based methods. Any damage to the system is theoretically reflected in changes to the vibration response, so exciting a system and monitoring its response allows quantification of damage, subject to changes in response being measurable and correlated to a specific damage type/location as necessary. In linear systems, there are many ways of implementing vibration-based structural health monitoring, principal among these being methods for analyzing changes in natural frequencies  and mode shapes [4,5]. Other, more exotic ways of monitoring changes involve examining flexibility , using time series analysis , or using wavelet analysis . Helpful reviews and details regarding modern computational algorithms for optimization and modeling are provided by Das et al.  and Gomes et al. . However, there are also many systems which are inherently nonlinear, or which become nonlinear once they have incurred damage. As an example, a beam which develops a breathing crack would, in many cases, no longer be able to be analyzed as a linear system. For these nonlinear systems, methods of linear analysis perform poorly and alternate analyses are currently limited , necessitating the development of additional methods.
Two linked issues critical to this endeavor are the choice of structural excitation and the method of processing the dynamic response. Currently, the most effective methods of nonlinear structural health monitoring excite (or interrogate) a structure using a chaotic signal and then measure and analyze the chaotic structural response via the machinery used for predicting chaotic time series to infer the presence and severity of damage. One of the earliest works in this area by Todd et al.  determined rules for how the Lyapunov exponents of the chaotic excitation should relate to the eigenvalues of the structure, and used a nonlinear time-prediction scheme called average local attractor variance ratio to detect damage. Researchers connected to Todd developed many other ways of detecting system changes involving fractal dimension , state space prediction , and nonlinear cross prediction error . These methods  were applied to mass-spring systems, beam systems, loosening bolts, and composite structures. A subset of these methods was also later revisited with hyperchaotic, instead of chaotic, excitation [17,18]. Other groups developed similar methods at about the same time, prominent among these being phase space warping [19–22] and sensitivity vector fields [23,24]. In general, this class of methods is effective and represents one of the few areas in which the search for practical applications of chaos theory succeeded. (Early in the study of chaos, investigators were anticipating many such applications  but few panned out.) The advantages of using chaotic interrogation, including broadband excitation, potential tuning of excitation to structural application, and sensitive dependence to initial conditions are well established.
What is less clear is whether methods used to analyze nonlinear dynamic response can be improved or extended. Nearly, all of the methods cited above adopt nonlinear time series prediction as an underlying algorithm in conjunction with a particular metric to assess system damage. While this type of analysis has proven effective at establishing the existence and quantifying the severity of damage, it has proven ineffective in determining the location or type of damage. Nonlinear time series prediction is poorly adapted to this task because a change in any system parameter influences the direction of exponential growth of the chaotic trajectories, making it difficult to discriminate between different parameter changes. As a result, only one parameter change can be reliably quantified at a time. Moreover, nonlinear time series prediction requires making many modeling choices, and it is usually unclear which will be best for a given application. The type of local model is one such choice (linear versus nonlinear neighborhoods, weighting schemes, etc.), and each model type typically also requires tuning additional parameters related to neighborhood size, time horizon, and other factors. These are often set using a trial-and-error process, which is why many of the studies cited above only make predictions one time-step ahead. Both the complexity and limitations of nonlinear time series prediction make it clear that complementary methods may be beneficial in order to extract the maximum useful information from a nonlinear structural response.
For this reason, several methods attempting to directly measure deformation of the system attractor have recently been developed. Instead of making trajectory predictions, these methods rely on measures of density or distribution to characterize the attractor and any changes it may undergo. The earliest attempt at this type of analysis appears to be by Liu and Chelidze  who evaluated system changes using the first moment of Poincare section data. Carroll and Byers [27,28] developed another method that compares the probability density of attractor data binned in histograms to infer system changes. Samadani et al.  also used density, but in a different way: estimating the density associated with a particular trajectory in a process termed phase space topology. Because, these methods are specifically designed for inferring system changes, they generally require less fine-tuning that their predictive counterparts. However, these methods are still limited to quantifying only one parameter change at a time.
In this paper, a new method is presented that uses chaotic excitation in combination with geometric, shape-based analysis of the boundary of a chaotic attractor to detect and quantify damage. Changes in the boundary of an attractor produce rich information about how parameters of the system are changing. The specific way this information is encoded is via boundary transformation vectors (BTVs). BTVs capture how parts of the boundary of an attractor must morph to align with the boundary of another, deformed attractor. Because this method works directly with attractor shape, it is better able to capture information about different types or locations of damage. BTVs will change in particular ways for one type/location of damage, while changing in other ways for second different type/location of damage. Additionally, damage level can be inferred using regions of BTVs that scale roughly proportionally to the damage. This combination represents a major benefit compared to past methods of damage detection using chaotic interrogation, whether in the nonlinear time series prediction or geometric deformation classes.
The remainder of the paper presents work showing how BTVs in conjunction with chaotic interrogation can be used to infer system damage. The initial sections summarize key ideas for performing chaotic interrogation and provide an algorithm for computing BTVs. Additional sections exhibit examples of BTVs being used to infer damage in simulated mass-spring-damper and cantilever beam systems. A summary along with directions for future investigation close the paper.
2.1 Chaotic Interrogation.
Chaotic interrogation is the use of a chaotic signal to excite a structure. The structure serves as a filter for the signal, altering it based on the properties of the system and any system variations, such as structural damage. By measuring the response of the structure, a second chaotic signal that encodes structural changes can be captured and used to infer damage.
Here S is used to indicate these Lyapunov exponents belong to the structure.
Here E is used to indicate these Lyapunov exponents belong to the excitation.
This necessary overlap of the structural and excitation Lyapunov spectra in forming the combined spectrum of the system is illustrated in Fig. 1. If becomes increasingly negative and less than structural exponent then, the modal participation of the structure increases as the ith mode becomes involved in the dynamics. An excitation time-scaling parameter δ can be used to control the spread of the excitation Lyapunov exponents and hence the spectral overlap.
Obviously, if becomes larger than then, the dimension of the attractor will increase, which will result in more complicated attractor that may be more difficult to analyze.
which is equivalent to creating a plane where . In this work, Poincare sections are constructed using the successive maxima formulation for all analyses to provide a standard way of examining the attractor.
2.2 Boundary Transformation Vectors.
Boundary transformation vectors describe the deformation that the boundary of one attractor would need to undergo in order to match the boundary of a second attractor. To encode this information, the shapes of two attractor boundaries are compared using a geometrical, shape-based transformation developed by Belongie et al.  called the shape context that has been adapted to attractor boundary matching . By comparing cataloged BTVs for known system damage to BTVs for unknown damage, unknown damage can be estimated. An overview of this idea is presented in Fig. 2.
The detailed process for comparing two attractors represented by Poincare sections, one section constructed using a signal associated with the undamaged structure (the reference section) and one section constructed using a signal associated with the damaged structure (the deformed section) is as follows:
Select regions of the attractor to compare: In many instances, Poincare sections constructed from experimental data may consist of disconnected segments. Some of these segments may be sensitive to damage (exhibiting significant deformation) while others may be relatively insensitive. In such situations, selecting particular segments or regions to analyze can be advantageous; selected regions are termed point clouds of interest (POI).
Find and smooth the boundary: The boundary of a POI can be captured using MATLAB's boundary.m function, which is based on work by Akkiraju et al. . To use this function, a value of shrink factor s must be chosen. The shrink factor controls compactness of the boundary and can be chosen based on the density and arrangement of Poincare section points. Because the tangent to each point on the boundary is also important to the algorithm, a quadratic polynomial is fit to groups of points of size ngroup centered on each boundary point, to obtain both a smoothed point and a tangent at that point. This process has been shown to be relatively insensitive to the size of ngroup . Here ngroup = 5.
Sample the boundary: Points along the smoothed boundary points are sampled, being selected sequentially a minimum distance of 2R apart, where R is the radius of a “hard core.” This ensures a relatively uniform sampling of the boundary and also preserves ordering, which is helpful for boundary matching.
- (4)Calculate shape contexts: Sampled boundaries will be matched using the idea of shape context, which is calculated for each sampled point. The shape context is defined bywhere hi is a histogram of the relative locations of the q points making up the shape in comparison to a reference point pi. This histogram uses 12 angular bins and 5 log-radial bins in a log-polar space to encode location information.(11)
- (5)Perform correspondence matching: The cost of matching points belonging to similar shapes can be calculated using the formulawhere CS is a cost based on matching shape context histograms, CT is a cost based on matching boundary tangents, and is a weighting factor determining the relative contributions of the two costs. The shape context cost is given by(12)where and are the K-bin normalized histograms for points pi and qj which belong to the reference and deformed sections, respectively. The boundary tangent cost is given by(13)where θi and θj are angles between the horizontal and tangent vectors associated with points pi and qj. Given a set of costs Cij between all pairs of points pi on one section and qj on the other, one seeks to minimize the total matching cost(14)subject to the constraint that π is a permutation. The cost minimization problem to assign matching points is solved using MATLAB's function matchpairs.m based on work by Duff and Koster . If the cost matrix is not square, this function accepts a parameter ϵ specifying the cost of nonassignment or any row or column. For a normalized cost matrix in which the maximum value of any element is 1, ϵ can be chosen in the range 0.1–0.5 to produce many quality matches.(15)
- (6)Calculate transformation coefficients: A transformation T(x, y) between the reference section and the deformed section can be calculated using thin-plate splines via the formulawhere the thin-plate spline interpolants and are given by(16)with and . Conditions for determining the coefficients of and are as follows:(17)
and for .
and have square integrable second derivatives.
and minimize the bending energy of their respective curves.
Here, m is the number of correspondences, αS is a scale factor, is a scale-independent regularization parameter controlling the amount of smoothing, and If is the bending energy. With vi set to and in turn for and , respectively, this gives thin-plate spline interpolant coefficients for the regularized problem.
(7) Iterate steps 5 and 6: The initial correspondence matching will typically contain some errors that impact the transformation estimate. By iterating steps 5 and 6 a few times (here five times), these inaccuracies are largely eliminated.
Additional details about this process as well as a parametric study of the sensitivity of the algorithm to its parameters are available in Ref. .
In this section, two different simulated structures are analyzed using BTVs to illustrate their applicability to structural health monitoring.
3.1 A Mass-Spring-Damper System.
This spectrum is desirable because the positive Lyapunov exponent is small, allowing scale δ to be adjusted a significant amount before Eq. (8) is violated. Plots of Poincare sections for the undamaged system excited for different values of the scale δ are shown in Fig. 4. For one system mode participates, for two system modes participate, but for Eq. (8) is violated, leading to an integer increase in the dimension of the attractor.
3.2 A Cantilever Beam.
The forced Brusselator described by Eq. (21) was also used to excite this structure for similar reasons to those described for the mass-spring-damper system. The value of was used as the forcing and applied at the free end of the beam. In this case, for one system mode participates, for two system modes participate, but for Eq. (8) is violated leading to an integer increase in the dimension of the attractor. Plots of Poincare sections for the undamaged beam excited for different values of the scale parameter δ are shown in Fig. 6.
For the mass-spring-damper system, results were generated using two scale factors: δ = 4 corresponding to one structural mode participating and δ = 8 corresponding to two structural modes participating. For scale δ = 4, the Poincare section is broken into discrete segments, so a POI must be selected for capturing the boundary; here the upper left portion of the section has been selected as indicated in Fig. 4. For scale δ = 8, the Poincare section is continuous, so the entire section can be considered as the POI. Key parameters for generating BTVs were set as s = 0.8, , and for all boundaries related to this system.
Example BTVs are shown in Fig. 7 for scale δ = 4 for a change in stiffness k1 from 2.0 to 1.9, with the Poincare section generated from measurements of position x8 associated with mass m8. Here, the blue points belong to the reference Poincare section associated with the undamaged structure, the red points belong to the deformed Poincare section associated with the damaged structure, and the green lines are BTVs approximating the attractor deformation. Points that remain unmatched following the final iteration of correspondence matching are boxed. The overall number of points involved in each POI are also indicated on the figure.
Figure 8 shows the magnitude and direction of BTVs generated for a range of variation in stiffness k1 from 1.95 to 1.70 compared to the undamaged system. The independent axis is a parameterized location along the attractor boundary running from 0 to 1. Note that the BTVs generally show changes in magnitude proportional to the level of damage the system has incurred; here regions 0.3–0.4 and 0.6–0.7 would allow an estimate of damage to be made simply by examining the magnitude of the BTVs.
Similar figures can be generated for changes in system damping. Figure 9 shows example BTVs for scale δ = 4 for changes in the damping c1 from 0.075 to 0.150, with the Poincare section again generated from measurements of position x8. The difference between this figure and Fig. 7 is that the direction of the BTVs have changed due to the damage being incurred in a different system parameter.
In Fig. 10, the magnitude and direction of BTVs generated for a range of variation in damping c1 from 0.100 to 0.200 are plotted. Again, for large regions of the plot, changes in the magnitude of the BTVs are proportional to the damage level. Note also that the angle of these vectors is different from when stiffness was changed as shown in Fig. 8, allowing the type of parameter change to be differentiated using vector direction.
Boundary transformation vectors can also be generated using Poincare sections with additional structural modes participating. Figure 11 shows a plot equivalent to Fig. 9 but now with scale δ = 8 (associated with two structural modes participating).
Because higher values of δ spread out the excitation Lyapunov spectrum and allow more overlap with the structural spectrum, richer patterns can occur. Figure 12 shows the magnitude and direction of BTVs associated with various changes in k1 and c1 at scale δ = 8. Note that at this scale, the directional signature for a change in stiffness versus damping is very distinct. Also note that while regions of boundary relative location where BTV magnitude is proportional to damage are smaller, they are still present for both types of parameter change.
Determining the location of structural damage is also possible. Figure 13 shows the magnitude and direction of BTVs associated with changes in stiffness at two locations: k1 and k7. Here changes in the damage location give rise to distinct patterns in the directions of the BTVs. Taken together, Figs. 12 and 13 show that the BTVs are a rich descriptor of how the system attractor is deformed, allowing changes in both damage level and type/location to be inferred.
These results are not limited to the mass-spring-damper system. For the cantilever beam system, results were generated using two scale factors δ = 10 corresponding to one structural mode participating and δ = 30 corresponding to two structural modes participating. For both scale δ = 10 and scale δ = 30, the Poincare section is broken into discrete segments, so a POI must be selected for drawing the boundary; the upper left portion of the section at each scale has been selected as indicated in Fig. 6. Key parameters for generating BTVs were set as s = 0.8, , and for all boundaries related to this structure.
In Fig. 14, example BTVs for the cantilever beam are shown for scale δ = 10 for a loss of 20% of mass and stiffness at element 5 (80% of undamaged values). Here, the Poincare section is generated from a recording of the displacement at the free end of the beam z20 (element 20).
Figure 15 shows the magnitude and direction of BTVs generated for incremental losses of material at element 5 of the cantilever beam. Here again changes in the magnitude of the BTVs are proportional to the damage incurred for some regions of the plot. The gap in the data is due to the boundary around the POI not being completely tight and therefore including a section without points connecting the two ends of the upside-down v-shape shown in Fig. 14.
Figure 16 shows that if the scale is increased to δ = 30 and two structural modes participate, then both amount and location of damage can be quantified in this system as well. In this figure, simulated material loss occurs at either element 5 or element 15 of the beam, leading to changes in BTV magnitude proportional to level of damage but distinct BTV directions depending on where the damage occurs.
The results shown in the immediately previous section demonstrate that BTVs generated via chaotic interrogation are a useful tool for assessing damage in structural health monitoring applications. Although the BTVs are not exact, they are systematic, reproducible, and offer damage estimates that are easy to understand and interpret. Moreover, because deformation of the attractor boundary is a rich descriptor of changes in the structural dynamics, BTVs can provide information about both damage type/location and amount, something not previously possible for methods based on chaotic interrogation.
Examining Fig. 8, we provide a brief example of how the magnitude of the BTVs can be used to assess damage level according to the process outlined in Fig. 2. If data regarding parameter values and Poincare sections were available for the nominal system with and at least two of the variations, say and , then Poincare sections and BTVs could be generated and used to create the magnitude and direction plots for those cases. These plots would be examined carefully to determine which boundary relative locations have BTV magnitudes that scale proportionally to system damage. Here choosing boundary relative location 0.7 would make sense, although large regions of the plot would likely also be acceptable, and several locations might be sampled to reduce uncertainty in any estimate. For a new, unknown change in the structure represented by either (interpolation) or (extrapolation) the proportionality at these locations could be used to estimate the value of k1 directly from the BTV magnitudes. At boundary relative location 0.7 estimates based on the BTVs are (interpolation) and (extrapolation) to three digits of precision. These estimates match the actual levels of damage in the system. Similar damage estimates could be made using data from Figs. 10, 12, 13, 15, or 16. All of these magnitude and direction plots have at least some relative boundary locations where BTV magnitude changes are roughly proportional to damage level. In fact, using BTV magnitudes at proportional locations for different excitation scales, such as at scale δ = 4 (Fig. 8) and again at scale δ = 8 (Fig. 12 or 13) might be another way to sample and increase reliability of any estimate in addition to choosing relative boundary locations.
As previously mentioned, one of the main advantages of using BTVs is to gain information about the type/location of damage in addition to inferring the damage level. This is possible even for minimal overlap of the structural and excitation Lyapunov spectra as evidenced by the difference in the typical direction of the BTVs for changes in k1 (Fig. 8) or c1 (Fig. 10) in the mass-spring-damper system at scale δ = 4. However, when the excitation scale increases the spectral overlap, indications of type/location are often more pronounced. This is demonstrated in Fig. 12 where BTV direction differences for changes in k1 and c1 are very distinct for the mass-spring-damper system, in Fig. 13 where differences between k1 and k7 changes are also very clear for the mass-spring-damper system, and in Fig. 16 where material losses at elements 5 and 15 are clearly differentiated for the cantilever beam. These results again suggest that interrogating a structure at multiple scales may be a valuable strategy, not only for quantifying damage but also for identifying its type/location. Obviously, there are limits to the amount the scale factor can be increased before the Poincare section loses structure and becomes less descriptive, but the results here show that for typical systems there may be a significant range to explore.
There are other limitations to using BTVs for damage detection. For one, it is possible that that damage incurred by a structure may become large enough that the shape of the Poincare section fundamentally changes. This is the case for the cantilever beam system if the stiffness/mass at element 5 drops below about 0.75. This suggests that caution is warranted in any case where an extrapolation is made based on a previous library of known BTVs, as an unknown boundary shape change could lie just beyond any previously investigated cases. Another potential issue not addressed in this paper is noise. Noise will impact the results in any real system and may limit the accuracy of BTV-based estimates of damage. However, previous work examining BTVs generated from signals with additive noise for chaotic systems  suggests that using simple, wavelet-based denoising will allow signals with signal-to-noise ratios of 100 or greater to be used reliably.
6 Conclusions and Future Work
This article demonstrates that boundary transformation vector analysis allows both damage level and damage type/location to be determined in simulated structural health monitoring applications involving linear systems. The ability to estimate both the amount and location of damage using BTVs represents a significant advantage over past methods making use of chaotic interrogation. Because BTVs are rich descriptors of how the boundary of an attractor deforms, they can encapsulate multiple aspects of how the system dynamics change when a structure experiences damage.
Work is currently underway to quantify damage in experimental linear systems using BTVs. However, since there are many ways to conduct structural health monitoring in linear systems, the ultimate goal is develop BTVs as a robust method of estimating damage in nonlinear systems. To this end, future work will focus on simulations of nonlinear systems. Additionally, the boundary transformation vector method should be expanded and refined. In this paper, Poincare sections were constructed using successive maxima of a structural signal because this is convenient when considering experimental data. However, Poincare sections constructed this way are likely not optimal in terms of sensitivity to damage; reconstructing the attractor using time-delayed embedding and determining a method for choosing optimal sections might be one way to further improve the method. Furthermore, in this work a forced Brusselator system was selected as the system excitation because it had a desirable Lyapunov spectrum; however, other choices of chaotic interrogation might perform even better. Methods to design or choose chaotic excitation tailored to a specific structure would go a long way to improving this and other damage detection methods based on chaotic interrogation.