## Abstract

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.

## 1 Introduction

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 [3] and mode shapes [4,5]. Other, more exotic ways of monitoring changes involve examining flexibility [6], using time series analysis [7], or using wavelet analysis [8]. Helpful reviews and details regarding modern computational algorithms for optimization and modeling are provided by Das et al. [9] and Gomes et al. [10]. 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 [11], 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. [12] 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 [13], state space prediction [14], and nonlinear cross prediction error [15]. These methods [16] 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 [1922] 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 [25] 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 [26] 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. [29] 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 Methods

### 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.

For the dynamics of the chaotic interrogation and the modes of the structure to interact, allowing changes in any structural parameters to be reflected in the system response, the Lyapunov spectrum of the structure and the excitation need to overlap. The Lyapunov exponents of a system describe the exponential rates at which system perturbations increase or decrease in different state space directions as time advances: a positive exponent indicates expansion; a negative exponent indicates contraction. The Lyapunov exponents of a structure are the real parts of the structure's eigenvalues. For instance, a linear structure described in standard state space form
$x˙=Ax+Bu$
(1)
where $x$ is a vector of the states of the system, A is the state matrix, B is the input matrix, and $u$ contains the system excitation, will have Lyapunov exponents given by the real parts of the eigenvalues of A. For typical structures, the Lyapunov exponents are negative forming the spectrum
$0≥λ1S>λ2S>…>λNS$
(2)

Here S is used to indicate these Lyapunov exponents belong to the structure.

Any excitation will have its own Lyapunov spectrum. For an excitation to be chaotic, one Lyapunov exponent must be positive, and the remainder of the exponents will be either zero or negative. Thus, the Lyapunov spectrum for a chaotic excitation is
$λ1E>0≥λ2E…≥λME$
(3)

Here E is used to indicate these Lyapunov exponents belong to the excitation.

For the Lyapunov spectra to overlap, so that structural changes are reflected in the dynamic response, at minimum the most negative Lyapunov exponent of the excitation must be less than the largest Lyapunov exponent of the structure
$λME<λ1S$
(4)

This necessary overlap of the structural and excitation Lyapunov spectra in forming the combined spectrum of the system is illustrated in Fig. 1. If $λME$ becomes increasingly negative and less than structural exponent $λiS$ 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.

Fig. 1
Fig. 1
Close modal
In addition to this necessary requirement, it is also desirable, for the purposes of analysis, to have the attractor of the excited structure exist in a low-dimensional space. The dynamics of the excited structure will exist in a space whose dimension is described by the Kaplan–Yorke conjecture
$D=j+∑i=1jλi|λj+1|$
(5)
where D is the fractal dimension of the dynamics and j is the largest number of exponents of the combined spectra that can be added before the sum becomes negative. Here, the combined spectrum is given by
$λi=(λjS,j=1,2…N;λkE,k=1,2,…M)$
(6)
$such that λ1>λ2>…>λN+M$
(7)
Thus, if a chaotic signal is used to excite a structure, for the lowest-dimensional response, it is necessary that
$|λ1E|<|λ1S|$
(8)
Because, this ensures that the dimension of the resulting attractor has the form
$D=2+λ1E|λ1S|$
(9)

Obviously, if $λ1E$ becomes larger than $|λ1S|$ then, the dimension of the attractor will increase, which will result in more complicated attractor that may be more difficult to analyze.

One approach to examining attractors is via Poincare sections. These are 2D slices of the attractor that can be constructed by forming a plane in the state space and recording data points when the trajectories of the system cross the plane. Another, sometimes more convenient, way of creating Poincare sections is to record successive maxima of a measured signal and then populate the 2D plane using data points
$point(i)=[xlocalmax i+1,xlocalmax i]$
(10)

which is equivalent to creating a plane where $xi˙=0$. 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. [30] called the shape context that has been adapted to attractor boundary matching [31]. 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.

Fig. 2
Fig. 2
Close modal

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:

1. (1)

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).

2. (2)

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. [32]. 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 [31]. Here ngroup = 5.

3. (3)

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. (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 by
$hi(k)=#{q≠pi:(q−pi)∈bin(k)}$
(11)
where 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.
5. (5)
Perform correspondence matching: The cost of matching points belonging to similar shapes can be calculated using the formula
$Cij=(1−βS)CS(pi,qj)+βSCT(pi,qj)$
(12)
where CS is a cost based on matching shape context histograms, CT is a cost based on matching boundary tangents, and $βS=0.25$ is a weighting factor determining the relative contributions of the two costs. The shape context cost is given by
$CS(pi,qj)=12∑k=1K[hi(k)−hj(k)]2hi(k)+hj(k)$
(13)
where $hi(k)$ and $hj(k)$ 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
$CT(pi,qj)=12‖( cos θi sin θi)−( cos θj sin θj)‖$
(14)
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
$H(π)=∑iCS(pi,qπ(i))$
(15)
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 [33]. 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.
6. (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 formula
$T(x,y)=(fx(x,y),fy(x,y))$
(16)
where the thin-plate spline interpolants $fx(x,y)$ and $fy(x,y)$ are given by
$f(x,y)=a1+axx+ayy+∑i=1mwiU(‖(xi,yi)−(x,y)‖)$
(17)
with $U(r)=r2 log r2$ and $U(0)=0$. Conditions for determining the coefficients of $fx(xi,yi)$ and $fy(xi,yi)$ are as follows:
• $fx(xi,yi)=xi′$ and $fy(xi,yi)=yi′$ for $i=1,2,…,m$.

• $fx(x,y)$ and $fy(x,y)$ have square integrable second derivatives.

• $fx(x,y)$ and $fy(x,y)$ minimize the bending energy of their respective curves.

Belongie et al. [30] provide additional detail on these conditions. In the presence of noise or other variability, it is often advantageous to relax the exact transformation T(x, y) by minimizing
$H|f|=∑i=1m(vi−f(xi,yi))2+αS2λ0If$
(18)

Here, m is the number of correspondences, αS is a scale factor, $λ0=0.5$ is a scale-independent regularization parameter controlling the amount of smoothing, and If is the bending energy. With vi set to $xi′$ and $yi′$ in turn for $fx(x,y)$ and $fy(x,y)$, 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. [31].

## 3 Applications

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.

Mass-spring-damper systems are a common testbed for damage detection methods. Here, BTV analysis is demonstrated using an eight degree-of-freedom system shown in Fig. 3 and described by Eq. (1) with A and B matrices given by
$A=[0IM−1KM−1C]; B=[0M−1]$
(19)
where M, C, and K are the structure's 8 × 8 banded mass, damping, and stiffness matrices, respectively, and forcing F is applied to the final mass m8. Here, we adopt individual mass, damping, and stiffness values of $mi=0.01, ci=0.075$, and $ki=2.0$ for $i=1…8$ to be consistent with past work [12,17] using chaotic interrogation. Damage was simulated by decreasing the stiffness of various system elements ki in increments of 0.05 or increasing the damping of various system elements ci in increments of 0.025. The Lyapunov spectrum of the undamaged system is given by
$λiS=−0.128,−1.12,−2.98,−5.45,−8.19,−10.8,−13.0,−14.5$
(20)
Fig. 3
Fig. 3
Close modal
Any number of chaotic systems could be used to excite this system, but here the forced Brusselator system is chosen
$u1˙=δ(u12u2−(b−1)u1+a+A sin u3)$
(21)
$u2˙=δ(−u12u2+bu1)$
(22)
$u3˙=ω$
(23)
where δ is a parameter used to scale time and tune the excitation's Lyapunov exponents for desired overlap with the system and the other parameters are the most common values for the Brusselator $[a b A ω]=[0.4 1.2 0.05 0.8]$. The value of u1 was used as the forcing. Prior to tuning, the forced Brusselator's Lyapunov spectrum is given by
$λiE=0.0140,0,−0.262$
(24)

This spectrum is desirable because the positive Lyapunov exponent $λ1E=0.0140$ 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 $0.488<δ<4.29$ one system mode participates, for $4.29<δ<11.4$ two system modes participate, but for $δ>9.12$ Eq. (8) is violated, leading to an integer increase in the dimension of the attractor.

Fig. 4
Fig. 4
Close modal

### 3.2 A Cantilever Beam.

Another easily simulated system for testing damage detection methods is a cantilever beam. Here, a Euler–Bernoulli beam with one end fixed and one end free is discretized using 20 finite elements, as shown in Fig. 5. To simulate this system, Eq. (1) can again be used, but here the system is cast in modal form so as to have control over the damping ratio associated with each mode, and thus the A and B matrices are given by
$A=[0I[ωi2]2[ζiωi]]; B=[0M−1VT]$
(25)
where $[ωi2]$ and $2[ζiωi]$ are 40 × 40 diagonal matrices of the beam's natural frequencies squared and double its Lyapunov exponents, respectively. The term $M−1VT$ in the B matrix represents the inverse of the mass matrix of the system multiplied by the eigenvectors which transform the system into modal coordinates. Here values for the beam's density ρ = 1185 kg/m3 and modulus E =3.2 GPa were adopted to be consistent with acrylic material and dimensions of length L =0.4128 m, width w =0.1016 m, and thickness t =0.0056 m were chosen to be reasonable for a tabletop experiment. The modal damping ratio was for all modes with a value of $ζi=0.01$. Damage was simulated by reducing both the mass and stiffness of a given element by 5% consistent with loss of material at a given location along the beam. The Lyapunov spectrum of the beam without damage is given by
$λiS=−0.5481,−3.4351,−9.6185,−18.8493,−31.1625,…$
(26)
Fig. 5
Fig. 5
Close modal

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 $3u1$ was used as the forcing and applied at the free end of the beam. In this case, for $2.09<δ<13.1$ one system mode participates, for $13.1<δ<36.7$ two system modes participate, but for $δ>39.2$ 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.

Fig. 6
Fig. 6
Close modal

## 4 Results

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, $2R=0.02$, and $ϵ=0.2$ 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.

Fig. 7
Fig. 7
Close modal

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.

Fig. 8
Fig. 8
Close modal

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.

Fig. 9
Fig. 9
Close modal

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.

Fig. 10
Fig. 10
Close modal

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).

Fig. 11
Fig. 11
Close modal

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.

Fig. 12
Fig. 12
Close modal

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.

Fig. 13
Fig. 13
Close modal

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, $2R=0.0001$, and $ϵ=0.2$ 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).

Fig. 14
Fig. 14
Close modal

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.

Fig. 15
Fig. 15
Close modal

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.

Fig. 16
Fig. 16
Close modal

## 5 Discussion

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 $k1=2.00$ and at least two of the variations, say $k1=1.95$ and $k1=1.85$, 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 $k1=0.90$ (interpolation) or $k1=0.80$ (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 $k1=1.90$ (interpolation) and $k1=1.80$ (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 [31] 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.

## References

1.
Doebling
,
S. W.
,
Farrar
,
C. R.
,
Prime
,
M. B.
, and
Shevitz
,
D. W.
,
1996
, “
Damage Identification and Health Monitoring of Structural and Mechanical Systems From Changes in Their Vibration Characteristics: A Literature Review
,” Los Alamos National Laboratory, Los Alamos, NM, Technical Report No. LA-13070-MS.
2.
Doebling
,
S. W.
,
Farrar
,
C. R.
, and
Prime
,
M. B.
,
1998
, “
A Summary Review of Vibration-Based Damage Identification Methods
,”
Shock Vib. Dig.
,
30
(
2
), pp.
91
105
.10.1177/058310249803000201
3.
Salawu
,
O. S.
,
1997
, “
Detection of Structural Damage Through Changes in Frequency: A Review
,”
Eng. Struct.
,
19
(
9
), pp.
718
723
.10.1016/S0141-0296(96)00149-6
4.
Pandey
,
A. K.
,
Biswas
,
M.
, and
Samman
,
M. M.
,
1991
, “
Damage Detection From Changes in Curvature Mode Shapes
,”
J. Sound Vib.
,
145
(
2
), pp.
321
332
.10.1016/0022-460X(91)90595-B
5.
Maia
,
N. M. M.
,
Silva
,
J. M. M.
,
Almas
,
E. A. M.
, and
Sampaio
,
R. P. C.
,
2003
, “
Damage Detection in Structures: From Mode Shape to Frequency Response Function Methods
,”
Mech. Syst. Signal Process.
,
17
(
3
), pp.
489
498
.10.1006/mssp.2002.1506
6.
Pandey
,
A. K.
, and
Biswas
,
M.
,
1994
, “
Damage Detection in Structures Using Changes in Flexibility
,”
J. Sound Vib.
,
169
(
1
), pp.
3
17
.10.1006/jsvi.1994.1002
7.
Sohn
,
H.
, and
Farrar
,
C. R.
,
2001
, “
Damage Diagnosis Using Time Series Analysis of Vibration Signals
,”
Smart Mater. Struct.
,
10
(
3
), pp.
446
451
.10.1088/0964-1726/10/3/304
8.
Taha
,
M. M. R.
,
Noureldin
,
A.
,
Lucero
,
J. L.
, and
Baca
,
T. J.
,
2006
, “
Wavelet Transform for Structural Health Monitoring: A Compendium of Uses and Features
,”
Struct. Health Monit.
,
5
(
3
), pp.
267
295
.10.1177/1475921706067741
9.
Das
,
S.
,
Saha
,
P.
, and
Patro
,
S. K.
,
2016
, “
Vibration-Based Damage Techniques Used for Health Monitoring of Structures: A Review
,”
J. Civ. Struct. Health Monit.
,
6
(
3
), pp.
477
507
.10.1007/s13349-016-0168-5
10.
Gomes
,
G. F.
,
Mendez
,
Y. A. D.
,
Alexandrino
,
P. L.
,
da Cunha
,
S. S.
, and
Ancelotti
,
A. C.
,
2019
, “
A Review of Vibration Based Inverse Methods for Damage Detection and Identification in Mechanical Structures Using Optimization Algorithms and ANN
,”
Arch. Comput. Methods Eng.
,
26
(
4
), pp.
883
897
.10.1007/s11831-018-9273-4
11.
Worden
,
K.
,
Farrar
,
C. R.
,
Haywood
,
J.
, and
Todd
,
M.
,
2008
, “
A Review of Nonlinear Dynamics Applications to Structural Health Monitoring
,”
Struct. Control Health Monit.
,
15
(
4
), pp.
540
567
.10.1002/stc.215
12.
Todd
,
M. D.
,
Nichols
,
J. M.
,
Pecora
,
L. M.
, and
Virgin
,
L. N.
,
2001
, “
Vibration-Based Damage Assessment Utilizing State Space Geometry Changes: Local Attractor Variance Ratio
,”
Smart Mater. Struct.
,
10
(
5
), pp.
1000
1008
.10.1088/0964-1726/10/5/316
13.
Nichols
,
J. M.
,
Virgin
,
L. N.
,
Todd
,
M. D.
, and
Nichols
,
J. D.
,
2003
, “
On the Use of Attractor Dimension as a Feature in Structural Health Monitoring
,”
Mech. Syst. Signal Process.
,
17
(
6
), pp.
1305
1320
.10.1006/mssp.2002.1521
14.
Nichols
,
J. M.
,
Todd
,
M. D.
, and
Wait
,
J. R.
,
2003
, “
Using State Space Predictive Modeling With Chaotic Interrogation in Detecting Joint Preload Loss in a Frame Structure Experiment
,”
Smart Mater. Struct.
,
12
(
4
), pp.
580
601
.10.1088/0964-1726/12/4/310
15.
Todd
,
M. D.
,
Erickson
,
K.
,
Chang
,
L.
,
Lee
,
K.
, and
Nichols
,
J. M.
,
2004
, “
Using Chaotic Interrogation and Attractor Nonlinear Cross-Prediction Error to Detect Fastener Preload Loss in an Aluminum Frame
,”
Chaos
,
14
(
2
), pp.
387
399
.10.1063/1.1688091
16.
Olson
,
C. C.
,
Overbey
,
L. A.
, and
Todd
,
M. D.
,
2005
, “
Sensitivity and Computational Comparison of State-Space Methods for Structural Health Monitoring
,”
Proc. SPIE
5768
, pp. 241–252.10.1117/12.598894
17.
Torkamani
,
S.
,
Butcher
,
E. A.
,
Todd
,
M. D.
, and
Park
,
G.
,
2011
, “
Detection of System Changes Due to Damage Using a Tuned Hyperchaotic Probe
,”
Smart Mater. Struct.
,
20
(
2
), p.
025006
.10.1088/0964-1726/20/2/025006
18.
Torkamani
,
S.
,
Butcher
,
E. A.
,
Todd
,
M. D.
, and
Park
,
G.
,
2012
, “
Hyperchaotic Probe for Damage Identification Using Nonlinear Prediction Error
,”
Mech. Syst. Signal Process.
,
29
, pp.
457
473
.10.1016/j.ymssp.2011.12.019
19.
Chelidze
,
D.
,
Cusumano
,
J. P.
, and
Chatterjee
,
A.
,
2002
, “
A Dynamical Systems Approach to Damage Evolution Tracking, Part 1: Description and Experimental Application
,”
ASME J. Vib. Acoust.
,
124
(
2
), pp.
250
257
.10.1115/1.1456908
20.
Cusumano
,
J. P.
,
Chelidze
,
D.
, and
Chatterjee
,
A.
,
2002
, “
A Dynamical Systems Approach to Damage Evolution Tracking, Part 2: Model-Based Validation and Physical Interpretation
,”
ASME J. Vib. Acoust.
,
124
(
2
), pp.
258
264
.10.1115/1.1456907
21.
Chelidze
,
D.
, and
Cusumano
,
J. P.
,
2004
, “
A Dynamical Systems Approach to Failure Prognosis
,”
ASME J. Vib. Acoust.
,
126
(
1
), pp.
2
8
.10.1115/1.1640638
22.
Chelidze
,
D.
, and
Cusumano
,
J. P.
,
2006
, “
Phase Space Warping: Non-Linear Time Series Analysis for Slowly Drifting Systems
,”
Philos. Trans. R. Soc., A
,
364
(
1846
), pp.
2495
2513
.10.1098/rsta.2006.1837
23.
Epureanu
,
B. I.
, and
Hashmi
,
A.
,
2006
, “
Parameter Reconstruction Based on Sensitivity Vector Fields
,”
ASME J. Vib. Acoust.
,
128
(
6
), pp.
732
740
.10.1115/1.2346692
24.
Sloboda
,
A. R.
, and
Epureanu
,
B. I.
,
2013
, “
Sensitivity Vector Fields in Time-Delay Coordinate Embeddings: Theory and Experiment
,”
Phys. Rev. E
,
87
(
2
), p.
022930
.10.1103/PhysRevE.87.022903
25.
Abarbanel
,
H. D. I.
,
1996
,
Analysis of Observed Chaotic Data
,
Springer-Verlag
,
New York
.
26.
Liu
,
M.
, and
Chelidze
,
D.
,
2006
, “
Identifying Damage Using Local Flow Variation Method
,”
Smart Mater. Struct.
,
15
(
6
), pp.
1830
1836
.10.1088/0964-1726/15/6/037
27.
Carroll
,
T. L.
,
2015
, “
Attractor Comparisons Based on Density
,”
Chaos
,
25
(
1
), p.
013111
.10.1063/1.4906342
28.
Carroll
,
T. L.
, and
Byers
,
J. M.
,
2016
, “
Grid-Based Partitioning for Comparing Attractors
,”
Phys. Rev. E
,
93
(
4
), p.
042206
.10.1103/PhysRevE.93.042206
29.
,
M.
,
Kwuimy
,
C. A. K.
, and
Nataraj
,
C.
,
2016
, “
Characterization of the Nonlinear Response of Defective Multi-DOF Oscillators Using the Method of Phase Space Topology (PST)
,”
Nonlinear Dyn.
,
86
(
3
), pp.
2023
2034
.10.1007/s11071-016-3012-x
30.
Belongie
,
S.
,
Malik
,
J.
, and
Puzicha
,
J.
,
2002
, “
Shape Matching and Object Recognition Using Shape Contexts
,”
IEEE Trans. Pattern Anal. Mach. Intell.
,
24
(
4
), pp.
509
522
.10.1109/34.993558
31.
Sloboda
,
A. R.
,
2021
, “
Boundary Transformation Representation of Attractor Shape Deformation
,”
Chaos
,
31
(
8
), p.
083133
.10.1063/5.0061029
32.
Akkiraju
,
N.
,
Edelsbrunner
,
H.
,
Facello
,
M.
,
Fu
,
P.
,
Mucke
,
E.
, and
Varela
,
C.
,
1995
, “
Alpha Shapes: Definition and Software
,”
Proceedings of the 1st International Computational Geometry Software Workshop
, Minneapolis, MN, Jan., p. 63–
66
.
33.
Duff
,
I. S.
, and
Koster
,
J.
,
2001
, “
On Algorithms for Permuting Large Entries to the Diagonal of a Sparse Matrix
,”
SIAM J. Matrix Anal. Appl.
,
22
(
4
), pp.
973
996
.10.1137/S0895479899358443