# Monitoring of the Looseness in Cargo Bolts under Random Excitation Based on Vibration Transmissibility

Authors: Guangming Dong; Jin Chen; Fagang Zhao
Journal: Shock and Vibration
 (2021)

---

## Abstract

Bolted joints are widely used in industrial applications and joint failure can cause a disastrous accident if loosening happens. Bolt loosening detection can be made by regular manual inspection or by using sensors based on different physical principles, such as acoustoelastic effect, piezoelectric active sensing, and electromechanical impedance. Compared with the above methods, vibration based bolt looseness monitoring using accelerometers is appealing for its economy and convenience for measurement. In this paper, cargo bolts looseness monitoring under random excitation is studied based on vibration transmissibility, which overcomes the drawback of commonly used vibration methods in finding local bolt looseness. Vibration transmissibility analysis only uses two vibration transducers to monitor bolt group looseness, where the vibration signal below the cargo bolts is used as the “input” and the other one above the cargo bolts is used as the “output.” There are 12 bolts in the cargo bolts studied in this paper, providing an essential clamping force to fix cargo during transportation. Six kinds of bolt group looseness with an increasing degree are simulated in the experiment. The experimental analysis shows that variation of the spectral moment can be used to monitor the global variation of the torque wrench exerted on the cargo bolts. The early stage of the bolt group looseness is that some one or two bolts begin to loose; however, the spectrum moment factor is insensitive to the local bolt looseness in the bolt group. To address this issue, the eigensystem realization algorithm (ERA) based on random input and output is utilized to find the subtle eigenvalue variation of the system matrix, which is neglected by the frequency transmissibility function. The experimental results show the effectiveness of the proposed method for detecting local bolt looseness.


---

## Body

## 1. Introduction

Bolted joints are widely used in mechanical and civil structures, such as machine tools, industrial robots, transport machines, bridges, and steel towers. A significant advantage of a bolted joint over other joint types, such as welded and riveted joints, is that they are capable of being disassembled. This feature, however, can cause joint failure leading to a disastrous accident if vibrational loosening happens, which is induced by flight/road loads or other environment uncertainties [1]. Hence, loosening detection of bolted joints is critical for structure safety, and it is highly desirable to achieve real-time online monitoring for bolts loosening detection under environmental and service conditions.The quality of a bolted connection is determined by bolt preload, which is usually controlled by the torque wrench method in engineering applications. However, it is difficult to achieve high precision of bolt preload because of the variation of friction between the bolt threads and the nut. Besides, the joints themselves may be located in inaccessible areas, which make regular inspection practically impossible.As a nondestructive testing technology, acoustoelastic effect based methods can be applied to determine the bolt preload stress by measuring the time difference of reflecting waves [2]. However, the time difference is only in the range of tens of nanoseconds, and the data acquisition process requires a high sampling rate for the data acquisition system. Besides, most references based on the acoustoelastic effect use oscilloscopes to capture the signals, which are not suitable for the in situ bolt health monitoring [3].The piezoelectric active sensing methods have been studied for bolt loosening detection. In the active sensing approach, one lead-zirconate-titanate (PZT) wafer is used as an actuator to generate an ultrasonic wave, and the other one is used as a sensor to record the wave that propagates through the bolted connection. The variation of the received wave energy [4–6], signal delay [7], spectral moment [8], and wavelet energy [9] have been experimentally investigated as an index for the bolting state. In the active sensing method, although the frequency for data acquisition can be lower than those used in the acoustoelastic effect based methods, it is still in the ultrasonic frequency band, which is quite a burden for data acquisition equipment in engineering applications.Besides the abovementioned active sensing methods, the electromechanical impedance (EMI) techniques have also been studied for bolt loosening detection [10–14]. EMI measures the electrical impedance of a PZT actuator-sensor patch bonded to a bolted structure. Due to electromechanical impedance coupling, the change in the structure properties such as stiffness and damping due to bolt loosening can be detected through the variation of the measured electrical impedance.The impedance method is attractive for local dynamics of the bolt joint with a relatively lower sensing frequency than the piezoelectric active sensing methods [3]. However, for effective monitoring of bolts group joints, a lot of PZT actuator-sensor patches are needed for bolt loosening detection. Besides, the impedance signal can also be altered by other ambient variations such as temperature changes, making it susceptible to false alarms [10].Vibration based bolt looseness monitoring using accelerometers is appealing for its economy and convenience for measurement. At present, much research about tool vibration based monitoring and fault diagnosis has been studied, which can be basically divided into two categories: (1) signal processing diagnostic method (2) and intelligent diagnostic method based on machine learning. Jiang and Yin generalized the recent advances in the multivariate statistical analysis based approaches for fault diagnosis and prognosis and developed a MATLAB toolbox for system fault detection (DB-KIT) [15]. Zhao et al. introduced an enhanced empirical wavelet transform for vibration based monitoring [16]; then they studied innovative machine learning methods for vibration based fault diagnosis [17] and prognosis [18]. The vibration based method is also effective for monitoring the connections with many bolted joints. The basic idea behind this technology is that any local damage somehow changes the dynamics of the structure and that the modal parameters identified from the structural response can reflect this change [19]. However, the measured vibration responses always reflect the global behavior of the whole structure. If there is another damage in the system that is under monitoring, or if the temperature and boundary conditions change, the bolt looseness monitoring based on vibration response will be inaccurate. To reduce the environmental effects, Cha et al. carried out a series of research on vision-based methods for bolt looseness detection, such as the Hough transform with support vector machines [20], the Viola–Jones algorithm with support vector machines [21], or convolutional neural networks [22, 23], and these studies can effectively reduce the impact of environmental changes on detection accuracy.Despite the limitations of the abovementioned vibration based techniques for bolt looseness monitoring, if there is a controlled input near the bolted structure, the vibration response can effectively reflect the bolt looseness. Razi et al. utilized an electric impact hammer to excite the bolted structure and analyzed the measured vibration using the empirical mode decomposition to establish an energy-based damage index [24]. With the impact hammer test, He and Zhu experimentally verified their study, which used the changes in natural frequencies to detect loosening of bolted connections in a full-size steel pipeline with bolted flanges [25]. Oregui et al. assessed the health condition of insulated rail joints in railway tracks by examining their dynamic response to impact excitation, where damage index is established by different frequency bands of vibration signals [26]. Wei et al. also utilized the exciting hammer to provide a stable excitation pulse signal and proposed a wavelet based damage index to monitor rail fastening status via a remote measurement system [27]. Kong et al. measured the responses of bolted joints, which were excited by impact, and employed the spectral features and the decision tree method to identify the bolt looseness level [28]. Using random excitation as a controlled input can be found in [29], where the authors studied statistical properties in time domain and signal power in different frequency bands for the measured vibration data, by which the bolt loosening can be detected.In many practices for online monitoring of structural health, the controlled input is difficult to be achieved. Thus, transmissibility-based structural health monitoring has attracted much attention by utilizing the output-to-output relationship [30, 31]. However, bolt looseness detection based on vibration transmissibility is rare in the literature, which is due to the fact that the responses on both sides of bolted connections may be uncorrelated at all. Fortunately, the cargo bolts in this study are another case: the motion of the whole structure is mainly in the vertical direction, and the responses below and above the bolted connection are obviously correlated.In this paper, the vibration signal below the cargo bolts is used as the “input,” and the vibration signal above the cargo bolts is used as the “output”; the concept of vibration transmissibility is used for bolt looseness monitoring. The frequency transmissibility function is calculated, and the derived spectrum moment factor is found to be related to the variation of torque level of the whole bolt group. However, the spectrum moment factor is insensitive to the local bolt looseness in the bolt group. To address this issue, the algorithm for the identification of a balanced state-space representation is utilized to find the subtle eigenvalue variation of the system matrix, which is neglected by the frequency transmissibility function. The first step of the proposed algorithm is the use of Markovsky’s algorithm [32] to estimate a finite part of the impulse response from input-output data under random environmental excitation during cargo transportation. Once the impulse response is computed, the second step is to obtain a balanced state-space model via the eigensystem realization algorithm (ERA) [33], by which the eigenvalues of the identified system matrix can be calculated as the bolt looseness indicator.The organization of this paper is as follows: the second section describes the studied cargo bolts and bolt looseness simulation; the third section shows the effectiveness of the frequency transmissibility function on bolt group looseness while also showing its insensitivity to local bolt looseness; and the fourth and fifth sections give the solution based on identification in the state-space domain.


## 2. The Studied Cargo Bolts and Looseness Simulation in Experiment

The cargo bolts are used in a clamping support structure, which is designed to fix cargo during transportation, as shown in Figure1. The enlarged view and cross section view of the structure are shown in Figures 2 and 3. The clamping support structure has 8 components as shown in Figure 3. The vertical braces (components 2 and 3) are fixed on the base (component 1). The vertical braces are also connected with two half-cylindrical clamps (components 4 and 5) by the cargo bolts (component 6 with 12 bolts, each side 6). When the cargo bolts are tightened, the simulated cargo (components 7 and 8, two half steel tubes) will be fixed by clamping force. The side view of the cargo bolts is also shown in Figure 2. Two vibration accelerometers below and above the cargo bolts are fixed on the structure.Figure 1
A clamping support structure used to fix cargo during transportation.Figure 2
Enlarged view of the clamping support structure with cargo bolts.Figure 3
Cross section view of the clamping support structure with cargo bolts.The studied cargo bolts in Figure2 have twelve bolts, evenly spaced along the plane connecting boundary. Under stochastic base vibration excitation, the early stage of the bolt group looseness is that some one or two bolts begin to loose, which increases dynamic loads on adjacent bolts, leading to the increased likelihood of bolt looseness. When more bolts loose, the previous loose bolts may become tight again due to the elasticity of the connecting boundary. However, the pretightening force of all the bolts in the bolt group will decrease. Thus, six damage cases of the bolt group looseness are simulated in the experiment: ① 2 bolts at one corner are loosened while the other 10 bolts are tightened; ② 4 bolts at the four corners are loosened while the other 8 bolts are tightened; ③ 4 bolts at corners are tightened, while the left 8 bolts are loosened; ④ the torque wrench levels of all the 12 bolts are decreased with a slight degree; ⑤ the torque wrench levels of all the 12 bolts are decreased with a medium degree; ⑥ the torque wrench levels of all the 12 bolts are decreased with a serious degree. The illustration of the first two damage cases is shown in Figure 4.Figure 4
Bolt looseness simulation for the first two cases.


## 3. The Proposed Algorithm for Cargo Bolts Looseness Monitoring

As discussed in Section2, the motion of the cargo bolts is mainly in the vertical direction; thus, the concept of vibration transmissibility is suitable for bolt looseness monitoring, where the vibration signal below the cargo bolts can be used as the “input” and the vibration signal above the cargo bolts are used as the “output.”In the following proposed algorithm, the frequency transmissibility function is calculated, and the derived spectrum moment factor is utilized to monitor the torque level variation of the whole bolt group. A brief theory description is shown in Section3.1 and corresponding experimental results are shown in Section 4.To further reflect the local bolt looseness in the bolt group, the algorithm for the identification of a balanced state-space representation is utilized to find the subtle eigenvalue variation of the system matrix. In the proposed state-space method, the first step is to estimate a finite part of the impulse response from input-output data under random environmental excitation during cargo transportation; the second step is to obtain a balanced state-space model via the eigensystem realization algorithm (ERA), by which the eigenvalues of the identified system matrix can be calculated as the local bolt looseness indicator. The theoretical background is introduced in Section3.2 and verification is shown in section 4.2.


### 3.1. Spectral Moment Factor

As shown in Figure2, two vibration accelerometers below and above the cargo bolts are fixed on the structure. In this paper, the vibration signal below the cargo bolts is used as the “input” and the vibration signal above the cargo bolts is used as the “output”; vibration transmissibility analysis is used for bolt looseness monitoring.At time instantk, let the sampled vibration signal below the cargo bolts be uk, and let the sampled signal above the cargo bolts be yk. The frequency transmissibility function (FTF) Sω can be calculated using the cross power spectrum Gyuω between uk and yk and using auto power spectrum Guuω of uk [31]:(1)Sω=GyuωGuuω.To quantify the cargo bolts looseness level, the spectral moment factor is defined as shown in equation (2), where only one input and one output are considered for simplicity. Spectral moments can describe the spectral distribution of a statistical signal, such as structural responses [34], EEG signals [35], and wave energy [36]. Spectral moments have been adopted for assessing the loosing state of a bolted structure by the piezoelectric active sensing method [8], where the requirements for additional energy consumption and high sampling frequency are needed. In this paper, spectral moments of vibration transmissibility function are calculated and used as an indicator of the torque level variation of the whole bolt group.Spectral moment with 0 and 1 order:(2)mi=∫0HωiSωdω,i=0,1.Spectral moment factor:(3)smf=m1m0.


### 3.2. Eigensystem Realization Algorithm (ERA) Based on Random Input and Output

#### 3.2.1. Estimation of the Impulse Response from Input-Output Data under Random Excitation

A state-space model for a dynamic system can be written as follows:(4)xk+1=Axk+Buk+wk,yk=Cxk+vk.The matrixA∈Rn×n is the state transition matrix that completely characterizes the dynamics of the system by its eigenvalues. The vectors uk∈Rm×1, yk∈Rp×1, and xk∈Rn×1 are the m inputs, p outputs, and n states, respectively, of the system at time instant k. To consider process noise and measurement noise of the model, noise vectors wk and vk are added into the model.Given the input vectoruk and the output vector yk, the state-space identification algorithms can be used to identify the system matrices A, B, and C, where the orthogonal or oblique projection of input and output data matrices are calculated through the QR decomposition, followed by Singular Value Decomposition (SVD) to determine the order, the observability matrix, and the state sequence. Then, the extraction of the state-space model is achieved through the solution of a least-squares problem. There are various kinds of state-space identification algorithms in the literature, such as N4SID, MOESP, and CVA, and Van Overschee incorporated them with different user-defined weighting matrices [37]. If only the estimation of the matrices A and C is concerned, from which modal parameters can then be extracted, the algorithms are generally equivalent [38].Markovsky proposed a computationally efficient algorithm for the computation of a finite part of the impulse responses directly from the input-output data [32]. Once the impulse response is computed, a balanced state model can be derived without returning to the original observed data.Suppose we haveN samples of input and output data: uk∈Rm×1, yk∈Rp×1, k∈0,…,N−1; the problem is how to compute the first L samples of impulse response matrices hk∈Rp×m, k∈0,…,L−1. The common sense is that N>L and L>nmax. nmax is an upper bound of the system order n.Similar to the definition of block Hankel matrices in [37], the following notation is used with f=f0,…,fN−1:(5)ℋΔf=f0f1…fN−Δf1f2…fN−Δ+1f2⋮⋮⋮⋮⋮fΔ−1fΔ…fN−1.Equation (5) defines the block Hankel matrix of f having Δ block rows. Then, through block Hankel matrices of input u and output y, matrices of Up, Uf, Yp, and Yf can be defined as follows:(6)ℋlmax+Lu=UpUf,ℋlmax+Ly=YpYf,where rowdimUp=rowdimYp=lmax, rowdimUf=rowdimYf=L, and row dim denotes the number of block rows of a matrix or vector. lmax is an upper bound of the system lag l.The impulse response sequence withL samples can be defined as a matrix with pL rows and m columns:(7)HpL×m=hT0,…,hTL−1T,where T denote the transpose of a matrix in equation (7). Clearly, HpL×m is the result of the system input as follows:(8)Uδ=Im0Lm−1×m,where I denotes the identity matrix with m rows and columns and 0 denotes the zero matrix with Lm−1 rows and m columns in equation (8).Under zero initial conditions, there exists a matrixG, such that(9)UpUfYpYfG=0mlmax×mUδ0plmax×mHpL×m.SinceG can be calculated through the first three submatrices in equation (9), the impulse sequences can be calculated as HpL×m=YfG. The calculation can be implemented via QR decomposition, and the details can be referred to in [32].


#### 3.2.2. ERA for System Matrix Estimation

Theoretically, the impulse sequences, or the Markov parameters, can be denoted using matricesA, B, and C defined in equation (4):(10)h0=CB,h1=CAB,…hk=CAkB,hk+1=CAk+1B,….ERA begins by forming thepα×mβ Hankel matrix ℋk composed of the impulse responses hk:(11)ℋk=hkhk+1…hk+βhk+1hk+2…hk+β+1hk+2⋮⋮⋮⋮⋮hk+αhk+α+1…hk+β+α.It can be derived that theℋk is connected with observability matrix Pα, the system matrix A , and controllability matrix Qβ:(12)ℋk=PαAkQβ,Pα=CCACA2⋮CAα,Qβ=BABA2B…AβB.The ERA process starts with the factorization of the Hankel matrixℋ0 by the use of SVD:(13)ℋ0=UΣVT,where the columns of matrices U and V are orthonormal and Σ is a rectangular matrix:(14)∑=Σr000,in which 0 denotes zero matrix with an appropriate dimension, and the matrix Σr is(15)Σr=diagσ1σ2…σiσi+1…σr,where σi is monotonically nonincreasing:(16)σ1≥σ2≥⋯≥σi≥σi+1≥⋯≥σr≥0.If we defineUr and Vr as the first r columns of the orthogonal matrices U and V, then we have(17)ℋ0=UrΣrVrT.Because of measurement noise, nonlinearity, and computer round-off error, the Hankel matrixℋ0 will usually be of full rank, which does not, in general, equal the true order of the system under test. By using SVD, noise and other interferences are removed through the truncation of nonzero small singular values.A balanced realization can be obtained by choosing the truncated observability matrix and controllability matrix as follows:(18)Pα=UrΣr1/2,Qβ=Σr1/2VrT.To compute the system matrixA, we can use the Hankel matrix ℋ1:(19)ℋ1=PαAQβ=UrΣr1/2AΣr1/2VrT.Then, the state transition matrixA can be calculated as(20)A=Σr−1/2UrTℋ1VrΣr−1/2.


#### 3.2.3. Physical Mode Determination through Stabilization Diagram and Clustering

The estimated system matrixA can be decomposed by eigendecomposition as A=ΨΛΨ−1, where Ψ is the eigenvector matrix and Λ=diagλdi is the diagonal eigenvalue matrix. The matrix Λ contains the n/2 discrete-time eigenvalues λdi and λdi∗ in complex conjugated pairs, which are directly related to the frequency and damping properties of the system. The discrete-time eigenvalues are first converted to continuous time eigenvalues λci as λci=lnλdi/dt , where dt is the time step of the digital data acquisition system. The undamped natural frequencies fni and damping ratios ζi can then be easily calculated from the conjugate pair of complex-valued eigenvalues: λci, λci∗=−2πζifni±j2πfni1−ζi2 .As the true system order is often unknown, a common practice is to calculate the modal parameters for increasing model ordersn. Besides, the weakly excited modes, which appear only at high model orders, can be detected. But if the model order n is higher than the true system order, the irrelevant noise will be modeled. These irrelevant noise modes can be visually removed in the stabilization diagram, in which the eigenfrequencies are plotted against the increasing model orders.The task can be very tedious in cases when there are closely coupled modes, heavily damped modes, and so forth. This creates a need for automating this selection process, thus avoiding user judgment by formulating a more reliable and intelligent process. A clustering algorithm is adopted to select the stable modes from the large set of parameters identified at various model orders. In this approach, one single set of modal estimates (a frequency and the associated damping ratio) is chosen among all the available parameters identified, and all sets that are closer to this point than a certain “radius” are clustered based on a distance measure between frequencies and damping estimates.The distance between modei and j is a Euclid distance that incorporates the eigenfrequency and damping ratio difference:(21)di,j=ωi−ωj2+ζi−ζj2.When certain clusters are selected, the variance of the clusters should be checked. When the variances of eigenfrequencies or damping ratios in each cluster are larger than a user-defined threshold, a conclusion of unsuccessful clustering will be made, and a new modal identification procedure should be restarted.


### 3.3. Flowchart of the Proposed Algorithm for Cargo Bolts Looseness Monitoring

Based on the theoretical description in Sections3.1 and 3.2, the flowchart of the cargo bolts looseness monitoring is shown in Figure 5. The proposed algorithm has four major steps:(1)
Cargo bolts are tightened for transportation. Because the cargo weight and bolt connection stiffness influence the calculated frequency transmissibility function, calibration is needed to calculate a baseline and threshold for a specific cargo with the tightened bolt connections. The algorithm for calibration is the same in the following monitoring procedure. The recalibration process for each transportation also helps to reduce the environmental impact, because the time for each transportation is relatively short.(2)
During cargo transportation, vibration signals below and above the bolted connection are measured as the “input” and “output” data for transmissibility analysis. The measurement is repeated for a constant duration. In this application, the duration is set as four seconds. The most time-consuming part is the estimation of the impulse response from input-output data, which is around 0.7 seconds in the experimental data processing. Other parts of the monitoring algorithm are less than 0.3 seconds. If the consuming times for consecutive data acquisition, display, and storage are included, the four-second duration is reasonable for online monitoring.(3)
During each four-second cycle, the spectral moment factor is calculated to monitor the torque level variation of the whole bolt group. A significant decrease of the spectral moment factor reflects a severe bolt loosening. When the factor is below the preset threshold, an emergency stop should be made for cargo bolts retightening.(4)
During each four-second cycle, the subtle eigenvalue variation of the system is calculated as the local bolt looseness indicator. When an abnormal frequency appears, while the spectral moment factor shows no noticeable change, a warning will arise in the monitoring system. In this case, the transportation goes on as usual until the next scheduled stop, where bolts connection state will be checked and retightened.Figure 5
Flowchart of the proposed algorithm for cargo bolts monitoring.


## 3.1. Spectral Moment Factor

As shown in Figure2, two vibration accelerometers below and above the cargo bolts are fixed on the structure. In this paper, the vibration signal below the cargo bolts is used as the “input” and the vibration signal above the cargo bolts is used as the “output”; vibration transmissibility analysis is used for bolt looseness monitoring.At time instantk, let the sampled vibration signal below the cargo bolts be uk, and let the sampled signal above the cargo bolts be yk. The frequency transmissibility function (FTF) Sω can be calculated using the cross power spectrum Gyuω between uk and yk and using auto power spectrum Guuω of uk [31]:(1)Sω=GyuωGuuω.To quantify the cargo bolts looseness level, the spectral moment factor is defined as shown in equation (2), where only one input and one output are considered for simplicity. Spectral moments can describe the spectral distribution of a statistical signal, such as structural responses [34], EEG signals [35], and wave energy [36]. Spectral moments have been adopted for assessing the loosing state of a bolted structure by the piezoelectric active sensing method [8], where the requirements for additional energy consumption and high sampling frequency are needed. In this paper, spectral moments of vibration transmissibility function are calculated and used as an indicator of the torque level variation of the whole bolt group.Spectral moment with 0 and 1 order:(2)mi=∫0HωiSωdω,i=0,1.Spectral moment factor:(3)smf=m1m0.


## 3.2. Eigensystem Realization Algorithm (ERA) Based on Random Input and Output

### 3.2.1. Estimation of the Impulse Response from Input-Output Data under Random Excitation

A state-space model for a dynamic system can be written as follows:(4)xk+1=Axk+Buk+wk,yk=Cxk+vk.The matrixA∈Rn×n is the state transition matrix that completely characterizes the dynamics of the system by its eigenvalues. The vectors uk∈Rm×1, yk∈Rp×1, and xk∈Rn×1 are the m inputs, p outputs, and n states, respectively, of the system at time instant k. To consider process noise and measurement noise of the model, noise vectors wk and vk are added into the model.Given the input vectoruk and the output vector yk, the state-space identification algorithms can be used to identify the system matrices A, B, and C, where the orthogonal or oblique projection of input and output data matrices are calculated through the QR decomposition, followed by Singular Value Decomposition (SVD) to determine the order, the observability matrix, and the state sequence. Then, the extraction of the state-space model is achieved through the solution of a least-squares problem. There are various kinds of state-space identification algorithms in the literature, such as N4SID, MOESP, and CVA, and Van Overschee incorporated them with different user-defined weighting matrices [37]. If only the estimation of the matrices A and C is concerned, from which modal parameters can then be extracted, the algorithms are generally equivalent [38].Markovsky proposed a computationally efficient algorithm for the computation of a finite part of the impulse responses directly from the input-output data [32]. Once the impulse response is computed, a balanced state model can be derived without returning to the original observed data.Suppose we haveN samples of input and output data: uk∈Rm×1, yk∈Rp×1, k∈0,…,N−1; the problem is how to compute the first L samples of impulse response matrices hk∈Rp×m, k∈0,…,L−1. The common sense is that N>L and L>nmax. nmax is an upper bound of the system order n.Similar to the definition of block Hankel matrices in [37], the following notation is used with f=f0,…,fN−1:(5)ℋΔf=f0f1…fN−Δf1f2…fN−Δ+1f2⋮⋮⋮⋮⋮fΔ−1fΔ…fN−1.Equation (5) defines the block Hankel matrix of f having Δ block rows. Then, through block Hankel matrices of input u and output y, matrices of Up, Uf, Yp, and Yf can be defined as follows:(6)ℋlmax+Lu=UpUf,ℋlmax+Ly=YpYf,where rowdimUp=rowdimYp=lmax, rowdimUf=rowdimYf=L, and row dim denotes the number of block rows of a matrix or vector. lmax is an upper bound of the system lag l.The impulse response sequence withL samples can be defined as a matrix with pL rows and m columns:(7)HpL×m=hT0,…,hTL−1T,where T denote the transpose of a matrix in equation (7). Clearly, HpL×m is the result of the system input as follows:(8)Uδ=Im0Lm−1×m,where I denotes the identity matrix with m rows and columns and 0 denotes the zero matrix with Lm−1 rows and m columns in equation (8).Under zero initial conditions, there exists a matrixG, such that(9)UpUfYpYfG=0mlmax×mUδ0plmax×mHpL×m.SinceG can be calculated through the first three submatrices in equation (9), the impulse sequences can be calculated as HpL×m=YfG. The calculation can be implemented via QR decomposition, and the details can be referred to in [32].


### 3.2.2. ERA for System Matrix Estimation

Theoretically, the impulse sequences, or the Markov parameters, can be denoted using matricesA, B, and C defined in equation (4):(10)h0=CB,h1=CAB,…hk=CAkB,hk+1=CAk+1B,….ERA begins by forming thepα×mβ Hankel matrix ℋk composed of the impulse responses hk:(11)ℋk=hkhk+1…hk+βhk+1hk+2…hk+β+1hk+2⋮⋮⋮⋮⋮hk+αhk+α+1…hk+β+α.It can be derived that theℋk is connected with observability matrix Pα, the system matrix A , and controllability matrix Qβ:(12)ℋk=PαAkQβ,Pα=CCACA2⋮CAα,Qβ=BABA2B…AβB.The ERA process starts with the factorization of the Hankel matrixℋ0 by the use of SVD:(13)ℋ0=UΣVT,where the columns of matrices U and V are orthonormal and Σ is a rectangular matrix:(14)∑=Σr000,in which 0 denotes zero matrix with an appropriate dimension, and the matrix Σr is(15)Σr=diagσ1σ2…σiσi+1…σr,where σi is monotonically nonincreasing:(16)σ1≥σ2≥⋯≥σi≥σi+1≥⋯≥σr≥0.If we defineUr and Vr as the first r columns of the orthogonal matrices U and V, then we have(17)ℋ0=UrΣrVrT.Because of measurement noise, nonlinearity, and computer round-off error, the Hankel matrixℋ0 will usually be of full rank, which does not, in general, equal the true order of the system under test. By using SVD, noise and other interferences are removed through the truncation of nonzero small singular values.A balanced realization can be obtained by choosing the truncated observability matrix and controllability matrix as follows:(18)Pα=UrΣr1/2,Qβ=Σr1/2VrT.To compute the system matrixA, we can use the Hankel matrix ℋ1:(19)ℋ1=PαAQβ=UrΣr1/2AΣr1/2VrT.Then, the state transition matrixA can be calculated as(20)A=Σr−1/2UrTℋ1VrΣr−1/2.


### 3.2.3. Physical Mode Determination through Stabilization Diagram and Clustering

The estimated system matrixA can be decomposed by eigendecomposition as A=ΨΛΨ−1, where Ψ is the eigenvector matrix and Λ=diagλdi is the diagonal eigenvalue matrix. The matrix Λ contains the n/2 discrete-time eigenvalues λdi and λdi∗ in complex conjugated pairs, which are directly related to the frequency and damping properties of the system. The discrete-time eigenvalues are first converted to continuous time eigenvalues λci as λci=lnλdi/dt , where dt is the time step of the digital data acquisition system. The undamped natural frequencies fni and damping ratios ζi can then be easily calculated from the conjugate pair of complex-valued eigenvalues: λci, λci∗=−2πζifni±j2πfni1−ζi2 .As the true system order is often unknown, a common practice is to calculate the modal parameters for increasing model ordersn. Besides, the weakly excited modes, which appear only at high model orders, can be detected. But if the model order n is higher than the true system order, the irrelevant noise will be modeled. These irrelevant noise modes can be visually removed in the stabilization diagram, in which the eigenfrequencies are plotted against the increasing model orders.The task can be very tedious in cases when there are closely coupled modes, heavily damped modes, and so forth. This creates a need for automating this selection process, thus avoiding user judgment by formulating a more reliable and intelligent process. A clustering algorithm is adopted to select the stable modes from the large set of parameters identified at various model orders. In this approach, one single set of modal estimates (a frequency and the associated damping ratio) is chosen among all the available parameters identified, and all sets that are closer to this point than a certain “radius” are clustered based on a distance measure between frequencies and damping estimates.The distance between modei and j is a Euclid distance that incorporates the eigenfrequency and damping ratio difference:(21)di,j=ωi−ωj2+ζi−ζj2.When certain clusters are selected, the variance of the clusters should be checked. When the variances of eigenfrequencies or damping ratios in each cluster are larger than a user-defined threshold, a conclusion of unsuccessful clustering will be made, and a new modal identification procedure should be restarted.


## 3.2.1. Estimation of the Impulse Response from Input-Output Data under Random Excitation

A state-space model for a dynamic system can be written as follows:(4)xk+1=Axk+Buk+wk,yk=Cxk+vk.The matrixA∈Rn×n is the state transition matrix that completely characterizes the dynamics of the system by its eigenvalues. The vectors uk∈Rm×1, yk∈Rp×1, and xk∈Rn×1 are the m inputs, p outputs, and n states, respectively, of the system at time instant k. To consider process noise and measurement noise of the model, noise vectors wk and vk are added into the model.Given the input vectoruk and the output vector yk, the state-space identification algorithms can be used to identify the system matrices A, B, and C, where the orthogonal or oblique projection of input and output data matrices are calculated through the QR decomposition, followed by Singular Value Decomposition (SVD) to determine the order, the observability matrix, and the state sequence. Then, the extraction of the state-space model is achieved through the solution of a least-squares problem. There are various kinds of state-space identification algorithms in the literature, such as N4SID, MOESP, and CVA, and Van Overschee incorporated them with different user-defined weighting matrices [37]. If only the estimation of the matrices A and C is concerned, from which modal parameters can then be extracted, the algorithms are generally equivalent [38].Markovsky proposed a computationally efficient algorithm for the computation of a finite part of the impulse responses directly from the input-output data [32]. Once the impulse response is computed, a balanced state model can be derived without returning to the original observed data.Suppose we haveN samples of input and output data: uk∈Rm×1, yk∈Rp×1, k∈0,…,N−1; the problem is how to compute the first L samples of impulse response matrices hk∈Rp×m, k∈0,…,L−1. The common sense is that N>L and L>nmax. nmax is an upper bound of the system order n.Similar to the definition of block Hankel matrices in [37], the following notation is used with f=f0,…,fN−1:(5)ℋΔf=f0f1…fN−Δf1f2…fN−Δ+1f2⋮⋮⋮⋮⋮fΔ−1fΔ…fN−1.Equation (5) defines the block Hankel matrix of f having Δ block rows. Then, through block Hankel matrices of input u and output y, matrices of Up, Uf, Yp, and Yf can be defined as follows:(6)ℋlmax+Lu=UpUf,ℋlmax+Ly=YpYf,where rowdimUp=rowdimYp=lmax, rowdimUf=rowdimYf=L, and row dim denotes the number of block rows of a matrix or vector. lmax is an upper bound of the system lag l.The impulse response sequence withL samples can be defined as a matrix with pL rows and m columns:(7)HpL×m=hT0,…,hTL−1T,where T denote the transpose of a matrix in equation (7). Clearly, HpL×m is the result of the system input as follows:(8)Uδ=Im0Lm−1×m,where I denotes the identity matrix with m rows and columns and 0 denotes the zero matrix with Lm−1 rows and m columns in equation (8).Under zero initial conditions, there exists a matrixG, such that(9)UpUfYpYfG=0mlmax×mUδ0plmax×mHpL×m.SinceG can be calculated through the first three submatrices in equation (9), the impulse sequences can be calculated as HpL×m=YfG. The calculation can be implemented via QR decomposition, and the details can be referred to in [32].


## 3.2.2. ERA for System Matrix Estimation

Theoretically, the impulse sequences, or the Markov parameters, can be denoted using matricesA, B, and C defined in equation (4):(10)h0=CB,h1=CAB,…hk=CAkB,hk+1=CAk+1B,….ERA begins by forming thepα×mβ Hankel matrix ℋk composed of the impulse responses hk:(11)ℋk=hkhk+1…hk+βhk+1hk+2…hk+β+1hk+2⋮⋮⋮⋮⋮hk+αhk+α+1…hk+β+α.It can be derived that theℋk is connected with observability matrix Pα, the system matrix A , and controllability matrix Qβ:(12)ℋk=PαAkQβ,Pα=CCACA2⋮CAα,Qβ=BABA2B…AβB.The ERA process starts with the factorization of the Hankel matrixℋ0 by the use of SVD:(13)ℋ0=UΣVT,where the columns of matrices U and V are orthonormal and Σ is a rectangular matrix:(14)∑=Σr000,in which 0 denotes zero matrix with an appropriate dimension, and the matrix Σr is(15)Σr=diagσ1σ2…σiσi+1…σr,where σi is monotonically nonincreasing:(16)σ1≥σ2≥⋯≥σi≥σi+1≥⋯≥σr≥0.If we defineUr and Vr as the first r columns of the orthogonal matrices U and V, then we have(17)ℋ0=UrΣrVrT.Because of measurement noise, nonlinearity, and computer round-off error, the Hankel matrixℋ0 will usually be of full rank, which does not, in general, equal the true order of the system under test. By using SVD, noise and other interferences are removed through the truncation of nonzero small singular values.A balanced realization can be obtained by choosing the truncated observability matrix and controllability matrix as follows:(18)Pα=UrΣr1/2,Qβ=Σr1/2VrT.To compute the system matrixA, we can use the Hankel matrix ℋ1:(19)ℋ1=PαAQβ=UrΣr1/2AΣr1/2VrT.Then, the state transition matrixA can be calculated as(20)A=Σr−1/2UrTℋ1VrΣr−1/2.


## 3.2.3. Physical Mode Determination through Stabilization Diagram and Clustering

The estimated system matrixA can be decomposed by eigendecomposition as A=ΨΛΨ−1, where Ψ is the eigenvector matrix and Λ=diagλdi is the diagonal eigenvalue matrix. The matrix Λ contains the n/2 discrete-time eigenvalues λdi and λdi∗ in complex conjugated pairs, which are directly related to the frequency and damping properties of the system. The discrete-time eigenvalues are first converted to continuous time eigenvalues λci as λci=lnλdi/dt , where dt is the time step of the digital data acquisition system. The undamped natural frequencies fni and damping ratios ζi can then be easily calculated from the conjugate pair of complex-valued eigenvalues: λci, λci∗=−2πζifni±j2πfni1−ζi2 .As the true system order is often unknown, a common practice is to calculate the modal parameters for increasing model ordersn. Besides, the weakly excited modes, which appear only at high model orders, can be detected. But if the model order n is higher than the true system order, the irrelevant noise will be modeled. These irrelevant noise modes can be visually removed in the stabilization diagram, in which the eigenfrequencies are plotted against the increasing model orders.The task can be very tedious in cases when there are closely coupled modes, heavily damped modes, and so forth. This creates a need for automating this selection process, thus avoiding user judgment by formulating a more reliable and intelligent process. A clustering algorithm is adopted to select the stable modes from the large set of parameters identified at various model orders. In this approach, one single set of modal estimates (a frequency and the associated damping ratio) is chosen among all the available parameters identified, and all sets that are closer to this point than a certain “radius” are clustered based on a distance measure between frequencies and damping estimates.The distance between modei and j is a Euclid distance that incorporates the eigenfrequency and damping ratio difference:(21)di,j=ωi−ωj2+ζi−ζj2.When certain clusters are selected, the variance of the clusters should be checked. When the variances of eigenfrequencies or damping ratios in each cluster are larger than a user-defined threshold, a conclusion of unsuccessful clustering will be made, and a new modal identification procedure should be restarted.


## 3.3. Flowchart of the Proposed Algorithm for Cargo Bolts Looseness Monitoring

Based on the theoretical description in Sections3.1 and 3.2, the flowchart of the cargo bolts looseness monitoring is shown in Figure 5. The proposed algorithm has four major steps:(1)
Cargo bolts are tightened for transportation. Because the cargo weight and bolt connection stiffness influence the calculated frequency transmissibility function, calibration is needed to calculate a baseline and threshold for a specific cargo with the tightened bolt connections. The algorithm for calibration is the same in the following monitoring procedure. The recalibration process for each transportation also helps to reduce the environmental impact, because the time for each transportation is relatively short.(2)
During cargo transportation, vibration signals below and above the bolted connection are measured as the “input” and “output” data for transmissibility analysis. The measurement is repeated for a constant duration. In this application, the duration is set as four seconds. The most time-consuming part is the estimation of the impulse response from input-output data, which is around 0.7 seconds in the experimental data processing. Other parts of the monitoring algorithm are less than 0.3 seconds. If the consuming times for consecutive data acquisition, display, and storage are included, the four-second duration is reasonable for online monitoring.(3)
During each four-second cycle, the spectral moment factor is calculated to monitor the torque level variation of the whole bolt group. A significant decrease of the spectral moment factor reflects a severe bolt loosening. When the factor is below the preset threshold, an emergency stop should be made for cargo bolts retightening.(4)
During each four-second cycle, the subtle eigenvalue variation of the system is calculated as the local bolt looseness indicator. When an abnormal frequency appears, while the spectral moment factor shows no noticeable change, a warning will arise in the monitoring system. In this case, the transportation goes on as usual until the next scheduled stop, where bolts connection state will be checked and retightened.Figure 5
Flowchart of the proposed algorithm for cargo bolts monitoring.


## 4. Experimental Results

### 4.1. Spectral Moment Factor as Looseness Indicator of the Cargo Bolts

With the whole structure under vertical random excitation in the experiment, the calculated FTFs for cargo bolts under two states are shown in Figures6 and 7. It can be seen in the figures that the spectral components are basically the same, while the peak positions decreased when bolt looseness happens. This is because of the decreased connection stiffness due to bolt looseness.Figure 6
FTF for cargo bolts in normal condition.Figure 7
FTF for cargo bolts in damage case ②.Variation of the defined spectral moment factor with the looseness level of cargo bolts is shown in Figure8. The looseness level 0 means the cargo bolts in normal condition. Looseness levels 1–6 correspond with the six damage cases defined in Section 2. Because the spectral components are mainly located in 0–1000 Hz, as shown in Figures 6 and 7, H is set as 1000 Hz when calculating the spectral moments in equation (2).Figure 8
Variation of the spectral moment factor with the bolt looseness level.It can be seen from Figure8 that the decreased value of spectral moment factor corresponds with the increase of looseness level for the cargo bolts. The vibration transmissibility analysis can be an effective tool for monitoring the cargo bolts looseness. However, the spectral moment factor is insensitive to local bolt looseness of damage case ①. To further increase the sensitivity of detecting local bolt looseness, vibration transmissibility analysis in the state-space domain is introduced and verified in the following sections. The experiment has been carried out three times, and the results have good repeatability.


### 4.2. Local Bolt Looseness Detection via Subtle Eigenvalue Variation of the System Matrix

#### 4.2.1. Simulation Demonstration

A simulation example is given at first to give an explanation of the subsequent experimental results. Figure9 shows a two-DOF system with base random excitation. The acceleration measured on the base plate is used as the system input, while the acceleration measured on the top mass element is used as the system output. The two-DOF system used in this simulation is tuned to have two close spaced natural frequencies with a relatively high damping ratio. The two natural frequencies are set as 1 Hz and 1.025 Hz, and the damping ratio is set as 0.016. Uncorrelated white noises wk and vk shown in equation (4) are added in the input and output signals to include uncertainty effects. Equation (4) describes a generalized state-space model with multi-input and multioutput. In this study, the input and output are all simplified to one channel.Figure 9
A two-DOF system with base random excitation.Figure10 shows the estimated impulse response of the two-DOF system from the two measured accelerations under random excitation. The sampling frequency is 25 Hz, and signal-noise-ratio (SNR) is set as 20 dB.Figure 10
Estimated impulse response from input-output data under random excitation for the two-DOF system.The calculated FTF using equation (1) is shown in Figure 11. Only one peak frequency with 1 Hz can be found, while the other natural frequency of 1.025 Hz is lost. The lack of frequency resolution is due to the effect of damping, which cannot be solved by simply increasing the acquired data length. The result using state-space domain identification is shown in Figure 12. The identified eigenfrequencies are confirmed to be a true physical mode because they are stable when the system order is increasing. The identified eigenfrequencies are 1 and 1.024 Hz, which are very close to their actual value.Figure 11
FTF of the two-DOF system.Figure 12
Variation of the identified frequencies with model order.When selecting the physical mode from the candidate eigenfrequencies, which are calculated from the system matrixA with a relatively high order, cluster analysis in the frequency-damping map is used. The distribution of identified frequency-damping is shown in Figure 13, where the noise effect is compared. When there is no noise, the identified mode is quite repeatable. When SNR is 20 dB, the identified results are very close to their actual value. As the SNR is decreasing, the identified results deviate from the actual value, but the error is still small. For example, the identified eigenfrequencies are 1.006 Hz and 1.035 Hz, respectively, and the error is less than 1%.Figure 13
The distribution of identified frequency-damping with model order from 4 to 22.A comparative study is made for different impulse responses estimation. In this paper, the impulse responseshk used in ERA are estimated using the state-space model as shown in Section 3.2.1, while the commonly used method for h˜k estimation is the inverse Fourier transform (IFFT) from the frequency transmissibility function Sω.The results usingh˜k are shown in Figures 14 and 15, which can be compared with Figures 12 and 13, where the length of impulse responses and parameters used in ERA is the same. Because there is variance in power spectra estimation for random signals using equation (1) [39], the estimated h˜k by IFFT of Sω has a larger error than hk using the state-space model. There are three closely spaced frequencies for different model orders, as shown in Figure 14. Cluster analysis is carried out in the identified frequency-damping map of Figure 15. They are identified as one cluster, and 1.008 Hz with damping ratio 0.045 is identified using equation (21). The identified result cannot actually separate the two closely spaced modes of 1 Hz and 1.025 Hz.Figure 14
Variation of the identified frequencies with model order when usingh˜k by IFFT of Sω.Figure 15
The distribution of identified frequency-damping with model order from 4 to 22, when usingh˜k by IFFT of Sω.


#### 4.2.2. Local Bolt Looseness Detection via Subtle Eigenvalue Variation of the System Matrix

As we can see from Section4.1, the spectral moment factor is insensitive to local bolt looseness of damage case ①. The 0–400 Hz components of the FTFs for cargo bolts under damage cases ① and ② are shown in Figure 16 and Figure 17. It can be seen from Figure 16 that the estimated base frequency using FTF is 202.5 Hz for damage case ①, while the base frequency of healthy cargo bolts is 203.5 Hz shown in Figure 6. The difference is quite a nuance and will be neglected in actual monitoring. When the bole looseness level is increased from case ① to case ②, the base frequency in the calculated FTF is decreased to 197 Hz. The results of Figures 16 and 17 can explain the result of the spectral moment factor in Figure 8, where the spectral moment factor for case ① does not changes while the spectral moment factor for case ② has a visible decrease.Figure 16
FTF (below 400 Hz) for cargo bolts in damage case ①.Figure 17
FTF (below 400 Hz) for cargo bolts in damage case ②.Analysis using the algorithm in Section3.2 is carried out in this section. Firstly, the estimated impulse response for the cargo bolts in damage case ① is shown in Figure 18. To avoid noise effect, impulse response with very small amplitudes is cut off. The sampling frequency is 3840 Hz.Figure 18
Estimated impulse response for cargo bolts under random excitation in damage case ①.The second step is to identify the system matrixA and calculate its eigenvalues using the estimated impulse response. From Figure 19 and Figure 20, the identified eigenfrequencies for the cargo bolts in damage case ① can be confirmed to be 196 Hz, 203 Hz, and 230 Hz. The identified base frequency is 196 Hz, which cannot be detected using FTF in Figure 16.Figure 19
Variation of the identified frequencies (below 400 Hz) with model order for damage case ①.Figure 20
The distribution of identified frequency-damping with model order for damage case ①.The results of Figure20 show that the identified base mode with 196 Hz is quite concentrated in the frequency-damping map, which shows relatively high energy compared with the scattered result of the third mode with 230 Hz. The emerge of 196 Hz is induced by asymmetric stiffness, which is caused by the local bolt looseness in the corner as shown in Figure 4. However, the base mode of 196 Hz is too close to the second mode of 203 Hz, which covers the base mode in the frequency domain using FTF in Figure 16. The results of Figure 19 and Figure 20 show that the analysis of vibration transmissibility in the state-space domain can discover the subtle eigenvalue variation of the system matrix.When the bolt looseness level is increased from case ① to case ②, the identification results are shown in Figures21 and 22. The base mode decreases from 203 Hz to 197 Hz due to the stiffness loss. As the damage case ② is four-bolt looseness at four corners, shown in Figure 4, the asymmetric effect does not exist. Thus, there are only two strong and distinct modes, and the frequency-damping is quite concentrated in Figure 22.Figure 21
Variation of the identified frequencies (below 400 Hz) with model order for damage case ②Figure 22
The distribution of identified frequency-damping with model order damage case ②.A comparative study is also made for case ① in this section for different impulse responses estimation, as conducted in the simulation Section4.2.1. The results using h˜k are shown in Figures 23 and 24, which can be compared with Figures 19 and 20, where the calculation parameters are the same. From the comparison of Figures 19 and 23, it can be seen that the estimated h˜k by IFFT of Sω has a larger error than hk using the state-space model. There are also three frequency lines in Figure 23, but the lines have more variance than Figure 19. Cluster analysis is carried out in the identified frequency-damping map of Figure 24, where it can be seen that the points of the first mode have a large dispersion area. The first two modes are identified as one cluster, and only two modes are identified with frequencies 199.8 Hz and 231.1 Hz.Figure 23
Variation of the identified frequencies (below 400 Hz) with model order when usingh˜k by IFFT of Sω, for damage case ①.Figure 24
The distribution of identified frequency-damping with model order when usingh˜k by IFFT of Sω, for damage case ①


## 4.1. Spectral Moment Factor as Looseness Indicator of the Cargo Bolts

With the whole structure under vertical random excitation in the experiment, the calculated FTFs for cargo bolts under two states are shown in Figures6 and 7. It can be seen in the figures that the spectral components are basically the same, while the peak positions decreased when bolt looseness happens. This is because of the decreased connection stiffness due to bolt looseness.Figure 6
FTF for cargo bolts in normal condition.Figure 7
FTF for cargo bolts in damage case ②.Variation of the defined spectral moment factor with the looseness level of cargo bolts is shown in Figure8. The looseness level 0 means the cargo bolts in normal condition. Looseness levels 1–6 correspond with the six damage cases defined in Section 2. Because the spectral components are mainly located in 0–1000 Hz, as shown in Figures 6 and 7, H is set as 1000 Hz when calculating the spectral moments in equation (2).Figure 8
Variation of the spectral moment factor with the bolt looseness level.It can be seen from Figure8 that the decreased value of spectral moment factor corresponds with the increase of looseness level for the cargo bolts. The vibration transmissibility analysis can be an effective tool for monitoring the cargo bolts looseness. However, the spectral moment factor is insensitive to local bolt looseness of damage case ①. To further increase the sensitivity of detecting local bolt looseness, vibration transmissibility analysis in the state-space domain is introduced and verified in the following sections. The experiment has been carried out three times, and the results have good repeatability.


## 4.2. Local Bolt Looseness Detection via Subtle Eigenvalue Variation of the System Matrix

### 4.2.1. Simulation Demonstration

A simulation example is given at first to give an explanation of the subsequent experimental results. Figure9 shows a two-DOF system with base random excitation. The acceleration measured on the base plate is used as the system input, while the acceleration measured on the top mass element is used as the system output. The two-DOF system used in this simulation is tuned to have two close spaced natural frequencies with a relatively high damping ratio. The two natural frequencies are set as 1 Hz and 1.025 Hz, and the damping ratio is set as 0.016. Uncorrelated white noises wk and vk shown in equation (4) are added in the input and output signals to include uncertainty effects. Equation (4) describes a generalized state-space model with multi-input and multioutput. In this study, the input and output are all simplified to one channel.Figure 9
A two-DOF system with base random excitation.Figure10 shows the estimated impulse response of the two-DOF system from the two measured accelerations under random excitation. The sampling frequency is 25 Hz, and signal-noise-ratio (SNR) is set as 20 dB.Figure 10
Estimated impulse response from input-output data under random excitation for the two-DOF system.The calculated FTF using equation (1) is shown in Figure 11. Only one peak frequency with 1 Hz can be found, while the other natural frequency of 1.025 Hz is lost. The lack of frequency resolution is due to the effect of damping, which cannot be solved by simply increasing the acquired data length. The result using state-space domain identification is shown in Figure 12. The identified eigenfrequencies are confirmed to be a true physical mode because they are stable when the system order is increasing. The identified eigenfrequencies are 1 and 1.024 Hz, which are very close to their actual value.Figure 11
FTF of the two-DOF system.Figure 12
Variation of the identified frequencies with model order.When selecting the physical mode from the candidate eigenfrequencies, which are calculated from the system matrixA with a relatively high order, cluster analysis in the frequency-damping map is used. The distribution of identified frequency-damping is shown in Figure 13, where the noise effect is compared. When there is no noise, the identified mode is quite repeatable. When SNR is 20 dB, the identified results are very close to their actual value. As the SNR is decreasing, the identified results deviate from the actual value, but the error is still small. For example, the identified eigenfrequencies are 1.006 Hz and 1.035 Hz, respectively, and the error is less than 1%.Figure 13
The distribution of identified frequency-damping with model order from 4 to 22.A comparative study is made for different impulse responses estimation. In this paper, the impulse responseshk used in ERA are estimated using the state-space model as shown in Section 3.2.1, while the commonly used method for h˜k estimation is the inverse Fourier transform (IFFT) from the frequency transmissibility function Sω.The results usingh˜k are shown in Figures 14 and 15, which can be compared with Figures 12 and 13, where the length of impulse responses and parameters used in ERA is the same. Because there is variance in power spectra estimation for random signals using equation (1) [39], the estimated h˜k by IFFT of Sω has a larger error than hk using the state-space model. There are three closely spaced frequencies for different model orders, as shown in Figure 14. Cluster analysis is carried out in the identified frequency-damping map of Figure 15. They are identified as one cluster, and 1.008 Hz with damping ratio 0.045 is identified using equation (21). The identified result cannot actually separate the two closely spaced modes of 1 Hz and 1.025 Hz.Figure 14
Variation of the identified frequencies with model order when usingh˜k by IFFT of Sω.Figure 15
The distribution of identified frequency-damping with model order from 4 to 22, when usingh˜k by IFFT of Sω.


### 4.2.2. Local Bolt Looseness Detection via Subtle Eigenvalue Variation of the System Matrix

As we can see from Section4.1, the spectral moment factor is insensitive to local bolt looseness of damage case ①. The 0–400 Hz components of the FTFs for cargo bolts under damage cases ① and ② are shown in Figure 16 and Figure 17. It can be seen from Figure 16 that the estimated base frequency using FTF is 202.5 Hz for damage case ①, while the base frequency of healthy cargo bolts is 203.5 Hz shown in Figure 6. The difference is quite a nuance and will be neglected in actual monitoring. When the bole looseness level is increased from case ① to case ②, the base frequency in the calculated FTF is decreased to 197 Hz. The results of Figures 16 and 17 can explain the result of the spectral moment factor in Figure 8, where the spectral moment factor for case ① does not changes while the spectral moment factor for case ② has a visible decrease.Figure 16
FTF (below 400 Hz) for cargo bolts in damage case ①.Figure 17
FTF (below 400 Hz) for cargo bolts in damage case ②.Analysis using the algorithm in Section3.2 is carried out in this section. Firstly, the estimated impulse response for the cargo bolts in damage case ① is shown in Figure 18. To avoid noise effect, impulse response with very small amplitudes is cut off. The sampling frequency is 3840 Hz.Figure 18
Estimated impulse response for cargo bolts under random excitation in damage case ①.The second step is to identify the system matrixA and calculate its eigenvalues using the estimated impulse response. From Figure 19 and Figure 20, the identified eigenfrequencies for the cargo bolts in damage case ① can be confirmed to be 196 Hz, 203 Hz, and 230 Hz. The identified base frequency is 196 Hz, which cannot be detected using FTF in Figure 16.Figure 19
Variation of the identified frequencies (below 400 Hz) with model order for damage case ①.Figure 20
The distribution of identified frequency-damping with model order for damage case ①.The results of Figure20 show that the identified base mode with 196 Hz is quite concentrated in the frequency-damping map, which shows relatively high energy compared with the scattered result of the third mode with 230 Hz. The emerge of 196 Hz is induced by asymmetric stiffness, which is caused by the local bolt looseness in the corner as shown in Figure 4. However, the base mode of 196 Hz is too close to the second mode of 203 Hz, which covers the base mode in the frequency domain using FTF in Figure 16. The results of Figure 19 and Figure 20 show that the analysis of vibration transmissibility in the state-space domain can discover the subtle eigenvalue variation of the system matrix.When the bolt looseness level is increased from case ① to case ②, the identification results are shown in Figures21 and 22. The base mode decreases from 203 Hz to 197 Hz due to the stiffness loss. As the damage case ② is four-bolt looseness at four corners, shown in Figure 4, the asymmetric effect does not exist. Thus, there are only two strong and distinct modes, and the frequency-damping is quite concentrated in Figure 22.Figure 21
Variation of the identified frequencies (below 400 Hz) with model order for damage case ②Figure 22
The distribution of identified frequency-damping with model order damage case ②.A comparative study is also made for case ① in this section for different impulse responses estimation, as conducted in the simulation Section4.2.1. The results using h˜k are shown in Figures 23 and 24, which can be compared with Figures 19 and 20, where the calculation parameters are the same. From the comparison of Figures 19 and 23, it can be seen that the estimated h˜k by IFFT of Sω has a larger error than hk using the state-space model. There are also three frequency lines in Figure 23, but the lines have more variance than Figure 19. Cluster analysis is carried out in the identified frequency-damping map of Figure 24, where it can be seen that the points of the first mode have a large dispersion area. The first two modes are identified as one cluster, and only two modes are identified with frequencies 199.8 Hz and 231.1 Hz.Figure 23
Variation of the identified frequencies (below 400 Hz) with model order when usingh˜k by IFFT of Sω, for damage case ①.Figure 24
The distribution of identified frequency-damping with model order when usingh˜k by IFFT of Sω, for damage case ①


## 4.2.1. Simulation Demonstration

A simulation example is given at first to give an explanation of the subsequent experimental results. Figure9 shows a two-DOF system with base random excitation. The acceleration measured on the base plate is used as the system input, while the acceleration measured on the top mass element is used as the system output. The two-DOF system used in this simulation is tuned to have two close spaced natural frequencies with a relatively high damping ratio. The two natural frequencies are set as 1 Hz and 1.025 Hz, and the damping ratio is set as 0.016. Uncorrelated white noises wk and vk shown in equation (4) are added in the input and output signals to include uncertainty effects. Equation (4) describes a generalized state-space model with multi-input and multioutput. In this study, the input and output are all simplified to one channel.Figure 9
A two-DOF system with base random excitation.Figure10 shows the estimated impulse response of the two-DOF system from the two measured accelerations under random excitation. The sampling frequency is 25 Hz, and signal-noise-ratio (SNR) is set as 20 dB.Figure 10
Estimated impulse response from input-output data under random excitation for the two-DOF system.The calculated FTF using equation (1) is shown in Figure 11. Only one peak frequency with 1 Hz can be found, while the other natural frequency of 1.025 Hz is lost. The lack of frequency resolution is due to the effect of damping, which cannot be solved by simply increasing the acquired data length. The result using state-space domain identification is shown in Figure 12. The identified eigenfrequencies are confirmed to be a true physical mode because they are stable when the system order is increasing. The identified eigenfrequencies are 1 and 1.024 Hz, which are very close to their actual value.Figure 11
FTF of the two-DOF system.Figure 12
Variation of the identified frequencies with model order.When selecting the physical mode from the candidate eigenfrequencies, which are calculated from the system matrixA with a relatively high order, cluster analysis in the frequency-damping map is used. The distribution of identified frequency-damping is shown in Figure 13, where the noise effect is compared. When there is no noise, the identified mode is quite repeatable. When SNR is 20 dB, the identified results are very close to their actual value. As the SNR is decreasing, the identified results deviate from the actual value, but the error is still small. For example, the identified eigenfrequencies are 1.006 Hz and 1.035 Hz, respectively, and the error is less than 1%.Figure 13
The distribution of identified frequency-damping with model order from 4 to 22.A comparative study is made for different impulse responses estimation. In this paper, the impulse responseshk used in ERA are estimated using the state-space model as shown in Section 3.2.1, while the commonly used method for h˜k estimation is the inverse Fourier transform (IFFT) from the frequency transmissibility function Sω.The results usingh˜k are shown in Figures 14 and 15, which can be compared with Figures 12 and 13, where the length of impulse responses and parameters used in ERA is the same. Because there is variance in power spectra estimation for random signals using equation (1) [39], the estimated h˜k by IFFT of Sω has a larger error than hk using the state-space model. There are three closely spaced frequencies for different model orders, as shown in Figure 14. Cluster analysis is carried out in the identified frequency-damping map of Figure 15. They are identified as one cluster, and 1.008 Hz with damping ratio 0.045 is identified using equation (21). The identified result cannot actually separate the two closely spaced modes of 1 Hz and 1.025 Hz.Figure 14
Variation of the identified frequencies with model order when usingh˜k by IFFT of Sω.Figure 15
The distribution of identified frequency-damping with model order from 4 to 22, when usingh˜k by IFFT of Sω.


## 4.2.2. Local Bolt Looseness Detection via Subtle Eigenvalue Variation of the System Matrix

As we can see from Section4.1, the spectral moment factor is insensitive to local bolt looseness of damage case ①. The 0–400 Hz components of the FTFs for cargo bolts under damage cases ① and ② are shown in Figure 16 and Figure 17. It can be seen from Figure 16 that the estimated base frequency using FTF is 202.5 Hz for damage case ①, while the base frequency of healthy cargo bolts is 203.5 Hz shown in Figure 6. The difference is quite a nuance and will be neglected in actual monitoring. When the bole looseness level is increased from case ① to case ②, the base frequency in the calculated FTF is decreased to 197 Hz. The results of Figures 16 and 17 can explain the result of the spectral moment factor in Figure 8, where the spectral moment factor for case ① does not changes while the spectral moment factor for case ② has a visible decrease.Figure 16
FTF (below 400 Hz) for cargo bolts in damage case ①.Figure 17
FTF (below 400 Hz) for cargo bolts in damage case ②.Analysis using the algorithm in Section3.2 is carried out in this section. Firstly, the estimated impulse response for the cargo bolts in damage case ① is shown in Figure 18. To avoid noise effect, impulse response with very small amplitudes is cut off. The sampling frequency is 3840 Hz.Figure 18
Estimated impulse response for cargo bolts under random excitation in damage case ①.The second step is to identify the system matrixA and calculate its eigenvalues using the estimated impulse response. From Figure 19 and Figure 20, the identified eigenfrequencies for the cargo bolts in damage case ① can be confirmed to be 196 Hz, 203 Hz, and 230 Hz. The identified base frequency is 196 Hz, which cannot be detected using FTF in Figure 16.Figure 19
Variation of the identified frequencies (below 400 Hz) with model order for damage case ①.Figure 20
The distribution of identified frequency-damping with model order for damage case ①.The results of Figure20 show that the identified base mode with 196 Hz is quite concentrated in the frequency-damping map, which shows relatively high energy compared with the scattered result of the third mode with 230 Hz. The emerge of 196 Hz is induced by asymmetric stiffness, which is caused by the local bolt looseness in the corner as shown in Figure 4. However, the base mode of 196 Hz is too close to the second mode of 203 Hz, which covers the base mode in the frequency domain using FTF in Figure 16. The results of Figure 19 and Figure 20 show that the analysis of vibration transmissibility in the state-space domain can discover the subtle eigenvalue variation of the system matrix.When the bolt looseness level is increased from case ① to case ②, the identification results are shown in Figures21 and 22. The base mode decreases from 203 Hz to 197 Hz due to the stiffness loss. As the damage case ② is four-bolt looseness at four corners, shown in Figure 4, the asymmetric effect does not exist. Thus, there are only two strong and distinct modes, and the frequency-damping is quite concentrated in Figure 22.Figure 21
Variation of the identified frequencies (below 400 Hz) with model order for damage case ②Figure 22
The distribution of identified frequency-damping with model order damage case ②.A comparative study is also made for case ① in this section for different impulse responses estimation, as conducted in the simulation Section4.2.1. The results using h˜k are shown in Figures 23 and 24, which can be compared with Figures 19 and 20, where the calculation parameters are the same. From the comparison of Figures 19 and 23, it can be seen that the estimated h˜k by IFFT of Sω has a larger error than hk using the state-space model. There are also three frequency lines in Figure 23, but the lines have more variance than Figure 19. Cluster analysis is carried out in the identified frequency-damping map of Figure 24, where it can be seen that the points of the first mode have a large dispersion area. The first two modes are identified as one cluster, and only two modes are identified with frequencies 199.8 Hz and 231.1 Hz.Figure 23
Variation of the identified frequencies (below 400 Hz) with model order when usingh˜k by IFFT of Sω, for damage case ①.Figure 24
The distribution of identified frequency-damping with model order when usingh˜k by IFFT of Sω, for damage case ①


## 5. Conclusion

This paper gives a case for vibration based bolt group looseness monitoring. The method uses the concept of vibration transmissibility. The motion of the whole structure studied in this paper is mainly in a vertical direction; thus, the use of the responses below and above the bolted connection as the “input” and “output” data for transmissibility analysis is feasible for bolt looseness detection.The frequency transmissibility function is calculated, and the derived spectrum moment factor is found to be related to the variation of torque level of the whole bolt group, where the decreased value of spectral moment factor corresponds with the increase of looseness level for the cargo bolts. The early stage of the bolt group looseness is that some one or two bolts begin to loose; however, the spectrum moment factor is insensitive to the local bolt looseness in the bolt group. To address this issue, the algorithm for the identification of a balanced state-space representation is further utilized to find the subtle eigenvalue variation of the system matrix, which is neglected by the frequency transmissibility function. The experimental results can discover the subtle eigenvalue variation of the system matrix and show the effectiveness of the proposed method for detecting local bolt looseness.



---

*Source: 8841940-2021-05-07.xml*