Structural Damage Identification Based on Variable-Length Elements and an Improved Genetic Algorithm for Railway Bridges

Abstract

A new damage identification method is proposed to solve the problem of no correspondence between the element division form of the finite element model and the actual damage location. The three basic operators in the traditional genetic algorithm are improved, and the catastrophe and neighborhood search processes are introduced to enhance the local optimization ability of the algorithm. The train–rail–bridge coupling time-varying equation is established. Based on the dynamic response of the bridge under trainload, the damage index is constructed, and the corresponding objective function is given. Through a numerical example, the stability and convergence rate of the algorithm are statistically analyzed. The effects of noise, the number of measuring points, and train speed on the recognition results are discussed. The research results indicate that, even if the damage location is different from the element division form of the finite element model, this method can accurately locate the damage location, but it will affect the quantitative results to a certain extent. In addition, the convergence speed of this method is fast, and the computing efficiency is about 6.7 times that of the conventional one-time recognition method.

1. Introduction

With the rapid development of the national economy, civil engineering structures tend to be large-scale and complex. In long-term service, infrastructure such as bridges, tracks, and overhead contact systems will be damaged due to the aging of materials, corrosion, long-term load action, and the influence of fire, earthquake, and other natural disasters [1,2,3]. The damage shortens the service life of the engineering structure and threatens the safety of peoples’ lives and property. Finding the location of structural damage in time and evaluating the degree of damage can ensure the safe operation of the structure and reduce the loss of individuals’ property [4]. In addition, the evolution law of structural performance deterioration can be mastered, the existing design methods can be improved, and the development of structural engineering can be promoted. Therefore, the research on structural damage identification methods has become a hot topic in civil engineering [5].

The core problem of structural damage identification is to judge the occurrence, location, and degree of damage. Since the 1980s, more and more attention has been paid to the research of structural dynamic damage identification based on environmental vibration [6]. Furthermore, a variety of research methods based on structural vibration response have been developed, such as the modal identification method [7], signal processing method [8], model updating method [9], and intelligent optimization algorithm [10]. The principle behind these methods is that damage will change the structural characteristics, such as stiffness, mass, flexibility, and damping, and then change the dynamic response of the structure [11]. Evaluating structural damage based on structural vibration response has been widely developed in the past two decades. For instance, Yang et al. [12] put forward the core idea and process of bridge frequency indirect identification based on cross-bridge vehicle response for the first time. It has the advantages of solid mobility, high efficiency, and economy. It can be tested continuously and quickly without road closure and stopping operation, so it has been rapidly promoted and developed [13,14]. With the deepening of the research, OBrien et al. [15] proposed a damage identification method using the intrinsic mode function (IMF) corresponding to the speed component response of the passing vehicle. Even in bridge deck irregularity, the distribution of bridge damage is successfully identified. However, if the damage is minor, the bridge deck irregularity may cover up the damage characteristics, making the damage challenging to effectively identify. Subsequently, Zhang et al. [16] used a one-dimensional convolution neural network to extract the features of structural acceleration to solve this problem. This method can identify small changes in local structural stiffness and mass and is overwhelmingly sensitive to small changes in local stiffness and mass.

One of the crucial problems of structural damage detection based on vibration is constructing and extracting sensitive parameters from structural dynamic responses to identify initial damage [17]. Intelligent algorithms, such as artificial neural network (ANN), cuckoo search (CS), and genetic algorithm (GA), are effectively applied to a wide range of optimization problems [18]. These algorithms are based on global search technology and evolve to find the best solution [19]. Their core idea is to use the vibration response of structures to construct the mapping relationship between characteristic damage parameters and actual damage and achieve the purpose of structural damage identification through reasoning, analysis, and calculation. To this end, an intelligent optimization algorithm has also made significant progress in damage identification. Its development trend is mainly reflected in the combination of various methods and the improvement of algorithms. For instance, Mousavi et al. [20] used an ANN as a damage detection method to evaluate the performance of the Hilbert–Huang transform based on fully integrated empirical mode decomposition and adaptive noise technology. The experimental results show that the method can accurately identify and extract the damage features of the signal, which indicates that the CEEMDAN-HT-ANN model has excellent ability and robustness in locating damage location and dividing damage severity. Huang et al. [21,22] proposed a structural damage identification method based on an improved whale intelligent optimization algorithm. This method takes the modal parameters as the identification factor, and several numerical examples are used to show that the method can accurately identify the damage location and severity of different structures and can effectively improve the efficiency of damage identification. Tran-Ngoc et al. [23] proposed a new damage identification method that combines ANN and CS algorithms. In numerical examples, various single intelligent algorithms are used to identify structural damage, which shows that this method is superior to other single intelligent algorithms and has higher computational efficiency. Nick et al. [24] proposed a two-stage damage detection method. The damage index based on modal strain energy is used for damage detection and location in the first stage. In the second stage, an ANN quantifies the damage degree. The results show that this method accurately predicts steel girder bridges’ damage location and damage size.

The above research shows that structural damage identification is a process from preliminary location to damage quantification [25]. Generally speaking, the bridge is invariably subjectively divided for structural damage identification. Then, the element’s damage is preset, and an appropriate method is selected for identification and verification according to the relationship between the characteristic damage parameters and the actual damage [26,27,28]. These studies have solved the problem of structural damage identification to some extent. Nevertheless, there is a problem in these treatments; it is assumed that the location of the damage is precisely an element in the finite element model (FEM), and the location of the damage corresponds one-to-one to the element division form of the FEM. However, this is an ideal and exceptional condition. In practice, the bridge structure is damaged due to the train’s impact, a lack of material strength, and the interplay of material faults and environmental conditions. Because these factors are unpredictable, the location and magnitude of these damages are uncertain and cannot be established without a professional evaluation. As a result, the damage location is not guaranteed to correspond one-to-one to the element division form of the FEM, and it is necessary to consider the difference between them.

In order to solve the problem, a multi-times damage location method based on variable-length elements and an improved genetic algorithm (IGA) is proposed. The convergence speed and recognition effect of the algorithm are studied. The effects of noise, the number of measuring points, and train speed on the recognition rate are discussed. First, a slightly larger length element is used to divide the bridge. Then, the appropriate damage step size and range are selected for rough identification to lock the scope of the bridge damage location. Finally, half of the length is used to subdivide these possible damage locations, the range of variables is limited, and then, the IGA is used to identify them. This process is repeated until the desired accuracy is achieved.

2. Theoretical Background

2.1. Genetic Algorithm

A GA was first proposed by Holland [29], which combines the principle of biological evolution, optimization technology, and computer technology and solves the extreme value problem by simulating the natural selection and genetic evolution mechanism of organisms in nature. It is a self-organizing and adaptive artificial intelligence technology. Because the GA generally does not act directly on the problem’s solution space but uses some coding method of the solution to express the answer, choosing an appropriate coding method is of great significance to the performance and efficiency of the algorithm.

Integer coding is used to map the solution space of the original problem to the integer string space, perform the genetic operation on the integer string space, and finally, restore the result to its phenotype through the decoding process to evaluate the fitness. Compared with binary coding, the length of integer strings is significantly reduced, and the decoding process is much simpler. Considering that the value ranges of decision variables in many optimization problems are often inconsistent, the variables with different value ranges can be grouped and coded so that any solution is composed of multiple chromosomes containing extra information, such as:

$$X=\left\{\begin{array}{c}{a}_0,a_1,a_2,\dots ,a_m\\ b_0,b_1,b_2,\dots ,b_n\\ \cdots \cdots \end{array}\right.$$

(1)

where X is an individual in the population, a and b are subcode strings of different variables, and m and n are the dimensions of the corresponding subcode strings, respectively.

Unlike real coding, each gene of integer coding is a decimal integer. Genetic operations do not directly act on the solution space. All genetic operations are carried out on each subcode string, which is more in line with the rules of gene calculation. It is especially suitable for the problem of discrete variables. Another advantage is avoiding the infeasible solution resulting in binary coding [30]. Suppose that a discrete variable has nine values, and binary coding requires 4-bit codes to represent sixteen solutions, resulting in seven kinds of invalid subcode strings. In integer coding, the variable can be defined by 0 to 8, and the genetic operation can make the gene take a value in the range of 0 to 8 without producing an invalid solution.

In addition to choosing the proper coding technique, the forms of operators in the GA, such as selection, crossover, and mutation, are also vital. They are the key carriers to imitate reproduction, hybridization, and mutation in the process of natural selection and heredity and form the foundation of the tremendous searchability of GAs. The classic genetic algorithm’s local optimization ability is inadequate, and it is simple to cause issues such as premature convergence, random roaming, and degradation. In this paper, three basic operators are improved, and catastrophe and neighborhood search mechanisms are introduced.

2.1.1. Hybrid Parallel Selection

The selection operation is based on the survival of the fittest concept, picking high-quality children and discarding low-quality individuals based on fitness. Its primary goal is to avoid gene deletion while improving global convergence and computing efficiency.

In this paper, a hybrid parallel selection operator is adopted: Preserving the Optimal Strategy and Tournament Strategy. New individuals are constantly produced in the GA through genetic operations such as crossover and mutation. Although an increasing number of excellent individuals will be built with the evolution of the population, they may also destroy the most adaptable individuals in the current population because of the randomness of the selection, crossover, mutation, and other operations [31]. Therefore, Preserving the Optimal Strategy is adopted to ensure that the final result of the iteration is the individual with the highest fitness all the time. There will be no degradation, and it can make the algorithm converge quickly. Tournament Strategy is a selection method based on the relationship between individuals’ fitness. The basic idea is to randomly select several individuals at a time and select the individual with the highest fitness to inherit to the next generation population.

• Tournament Strategy

M individuals were randomly selected from the population to compare their fitness, and the best individual was chosen as the parent, repeated N times.

b.

Preserving the Optimal Strategy

The best parent from the current population is selected and copied directly to the next generation without crossover or mutation.

2.1.2. Hybrid Crossover

Crossover is the process of transferring genetic material from one individual to another to create superior children. It is a crucial element separating GAs from other evolutionary algorithms. Moreover, it is the primary mechanism for producing new individuals and plays a critical part in GAs.

The Hybrid Crossover operator employed in this study is derived from the One-Point Crossover operator in binary coding [32]. An intersection is randomly set in the individual coding string, and part of the chromosomes of the two paired individuals at the point are then exchanged. The operator performs head–head crossover and head–tail crossover, respectively, according to a certain probability of p_H and (1-p_H) to expand the search range and improve the ability of the algorithm to jump out of the local optimization. Suppose that the two fathers are A and B, taking the case of a single chromosome as an example. The schematic diagram of the new individuals ${\mathrm{A}}^{\prime }$ and ${\mathrm{B}}^{\prime }$ is shown in Figure 1.

2.1.3. Adjacent Mutation

Mutation is a supplementary approach for recovering an individual’s lost or underdeveloped genetic material to prevent children from converging prematurely in the process of producing the optimal solution. It is also an essential operation step to replace some gene values in the individual coding string with other gene values to form a new individual.

In this paper, the Adjacent Mutation operator is used. The positions of mutations are randomly determined in the integer code string. According to the preset mutation probability p_m, the gene values of these mutation positions are randomly added or subtracted by one. The gene value is only added by one if it is the minimum. Similarly, if it is the maximum, the value of the gene is only reduced by one.

$$b_k^{\prime }=\left\{\begin{array}{l}{b}_k+1,b_{k\min }<b_k<b_{k\max }\&r<0.5\\ b_k+1,b_k=b_{k\min }\\ b_k-1,b_{k\min }<b_k<b_{k\max }\&r\ge 0.5\\ b_k-1,b_k=b_{k\max }\end{array}\right.$$

(2)

where $b_k$ and $b_k^{\prime }$ are the values before and after a gene mutation, respectively, and the subscript max and min are the upper and lower limits of the gene value, respectively.

2.1.4. Catastrophe

Preserving the Optimal Strategy in the selection operator can make the algorithm converge quickly, but making the algorithm fall into local optimization is also easy. Hence, it is necessary to improve the algorithm’s global search performance further. If the initial population size increases, the algorithm’s performance will be improved, but it will also greatly increase the solving time [33]. Therefore, a revolutionary force is required to modify the centralized distribution of initial solutions to boost the GA’s global search performance.

The idea of catastrophe in biological evolution is introduced in this study. Significant changes in the external environment, such as ice ages, forest fires, earthquakes, and floods, are referred to as catastrophes. The destruction of the great majority of creatures, culminating in the extinction of many species, is commonly referred to as a catastrophe. Only species or individuals with high adaptation will be able to live and reproduce after the disaster. Species or individuals that have survived a “catastrophe” have a higher survival rate. The simple steps of the catastrophe operator are as follows:

• The number of catastrophes is initialized.
• After the crossover and mutation are carried out in the current evolution stage, the population fitness is evaluated. If the number of individuals with the same fitness exceeds the threshold, individual fitness is ranked.
• The random generation of m individuals replaces the individuals with lower fitness until the maximum number of catastrophes is reached.

2.1.5. Neighborhood Search

The population diversity decreases dramatically as the number of iterations grows. Genetic algorithms’ crossover and selection operators can no longer produce new individuals with greater vitality, and it is unlikely that the lost effective genes can be recovered solely through mutation. Therefore, the neighborhood search mechanism is introduced to further improve local searchability in this paper.

At the end of the iterative process, the optimal solution is added and subtracted by one, respectively, and the full permutation in mathematics is then carried out. All the combinations in the full permutation are evaluated for fitness and the best is selected. For instance, if [a b c d] is a set of optimal solutions, there are $3^4=81$ possibilities of combinations of these solutions.

2.1.6. Performance Evaluation

In order to evaluate the actual optimization performance of the IGA, two commonly used optimization algorithm test functions are introduced and compared with the traditional GA. The test functions are shown below [34].

$$f_1\left(x\right)={\displaystyle \sum _{i=1}^{20}{x}_i^2},x_i\in [-100,100] \& x_i\mathrm{is}\mathrm{an}\mathrm{integer}$$

$$f_2(x)={\displaystyle \sum _{i=1}^{30}\left|x_i\right|}+{\displaystyle \prod _{i=1}^{30}\left|x_i\right|}, x_i\in [-30,30]\&x_i\mathrm{is}\mathrm{an}\mathrm{integer}$$

The population size is set to 100, the maximum number of iterations is 500, and the initial population in the two algorithms is the same. The optimization calculation of $f_1\left(x\right)$ and $f_2\left(x\right)$ is carried out, respectively, and the objective function (OBJ) iteration curve in each offspring is obtained, as shown in Figure 2.

As shown in Figure 2, when searching for the optimal global solution, the convergence speed of the IGA is significantly higher than the GA, and the fitness value is less than the GA, which shows that the IGA is better than the GA in the optimization process.

2.2. Vertical Coupling Vibration Model of Train–Rail–Bridge

When a train is running on a bridge, the bridge will vibrate and deform due to the trainload, and the vibration deformation of the bridge will be fed back to the train, which will cause the vibration deformation of the running train. Therefore, in the problem of train–bridge coupling vibration caused by trains crossing, the dynamic response of the bridge and the train contains the modal or geometric parameter information of the bridge structure [12]. With the rapid development of vibration signal testing technology and the improvement of vehicle–rail–bridge coupling vibration theory, structural damage identification based on dynamic response has received widespread attention in recent years because of its advantages, such as strong mobility, high efficiency, and economy. It can detect continuously and quickly without road closure and parking operation, so it has a good application prospect [14].

The whole coupling system is divided into train–rail–bridge systems. The Euler–Bernoulli beam element with six degrees of freedom (DOFs) is used to model rail and simply supported beams [35], and the direct stiffness method is used to couple the rail and the bridge. Hertz’s contact theory connects the train and the rail system. Multi-rigid-body dynamics are used to model the train system with 10 DOFs [36]. When considering various coupling relations and irregularities, the structure is discretized, and the coupling elements of an arbitrary length combination, as shown in Figure 3, are obtained.

2.2.1. Train–Rail Vertical Coupling Relationship

The train–rail–bridge interaction relies on wheel–rail contact. It is not only the train and lower rail’s connection manner but also the two subsystems’ primary excitation source. A Hertz spring with the stiffness of k_H is used to imitate this contact in a mobile wheel–rail interaction element, as shown in Figure 3. According to the linearized Hertz elastic contact theory, the contact stiffness is calculated using the tangent slope method [37].

$$k_{\mathrm{H}}=\frac{\mathrm{d}[P(t)]}{\mathrm{d}[\delta \mathrm{z}(t)]}=\frac{3}{2G}{P}_0^{1/3}$$

(3)

where G is the wheel–rail contact constant (m/N^2/3), taking G = 3.86r^−0.115 × 10⁻⁸; R is the radius of the wheel rolling circle (m); P₀ is the static axle load of the train.

The standard elastic compression of the wheel–rail contact point is:

$$\delta z\left(t\right)=z_{\mathrm{w}}\left(t\right)-z_{\mathrm{r}}\left(x_{\mathrm{w}},t\right)-r\left(x_{\mathrm{w}}\right)$$

(4)

where $z_{\mathrm{w}}\left(t\right)$ is the vertical displacement of the wheel at t time, $z_{\mathrm{r}}\left(x_{\mathrm{w}},t\right)$ is the vertical displacement of the rail where the wheel acts on $x_{\mathrm{w}}$ at t time, $x_{\mathrm{w}}$ is the position coordinate of the wheel, and $r\left(x_{\mathrm{w}}\right)$ is the value of rail irregularity at the $x_{\mathrm{w}}$ of the wheel.

The vertical displacement of the rail at the contact point of the wheel can be obtained from the node displacement vector $z_{\mathrm{r}}^{\mathrm{e}}$ of the rail element and the corresponding cubic Hermite interpolation shape function, that is:

$$z_{\mathrm{r}}\left(x_{\mathrm{w}},t\right)=N_{x=x_{\mathrm{w}}}{z}_{\mathrm{r}}^{\mathrm{e}}$$

(5)

$$\delta z\left(t\right)=z_{\mathrm{w}}\left(t\right)-z_{\mathrm{r}}\left(x_{\mathrm{w}},t\right)-r\left(x_{\mathrm{w}}\right)$$

where N is the cubic Hermite interpolation vector of the displacement of the rail and bridge elements. The indirect contact force between the wheel and the rail can be expressed as [38]:

$$F_{\mathrm{in}}=-a_{\mathrm{in}}{k}_{\mathrm{H}}\delta z\left(t\right)$$

(6)

$$a_{\mathrm{in}}=\left\{\begin{array}{cc}1& ,\delta z\left(t\right)\ge 0\\ 0& ,\delta z\left(t\right)<0\end{array}\right.$$

2.2.2. Bridge Element Equation

The dynamic equation of the bridge element is:

$$M_{\mathrm{b}}^{\mathrm{e}}{\ddot{z}}_{\mathrm{b}}^{\mathrm{e}}+C_{\mathrm{b}}^{\mathrm{e}}{\dot{z}}_{\mathrm{b}}^{\mathrm{e}}+K_{\mathrm{b}}^{\mathrm{e}}{z}_{\mathrm{b}}^{\mathrm{e}}=P_{\mathrm{b}}^{\mathrm{e}}-{\displaystyle {\int }_0^{l_b}{N}_{\mathrm{b}}^{\mathrm{T}}{q}_{\mathrm{rb}}}\left(x,t\right)\mathrm{d}x$$

(7)

where $M_{\mathrm{b}}^{\mathrm{e}}$, $C_{\mathrm{b}}^{\mathrm{e}}$, and $K_{\mathrm{b}}^{\mathrm{e}}$ are bridge element mass matrix, damping matrix, and stiffness matrix, respectively; $z_{\mathrm{b}}^{\mathrm{e}}$ is the node displacement vector of bridge element; $P_{\mathrm{b}}^{\mathrm{e}}$ is the load vector of the element node, and $l_{\mathrm{b}}$ is the length of a bridge element.

The distributed force $q_{\mathrm{rb}}$ transmitted from the rail to the bridge is:

$$q_{\mathrm{rb}}\left(x,t\right)=-k_{\mathrm{rb}}\left(z_{\mathrm{pr}}-z_{\mathrm{b}}\right)-c_{\mathrm{rb}}\left({\dot{z}}_{\mathrm{pr}}-{\dot{z}}_{\mathrm{b}}\right)$$

(8)

$$\begin{gathered} z_{p\mathrm{r}}=N_{p\mathrm{r}}{z}_{p\mathrm{r}}^{\mathrm{e}},{\dot{z}}_{p\mathrm{r}}=N_{p\mathrm{r}}{\dot{z}}_{p\mathrm{r}}^{\mathrm{e}} \\ {z}_{p\mathrm{r}}=N_{p\mathrm{r}}{z}_{p\mathrm{r}}^{\mathrm{e}},{\dot{z}}_{p\mathrm{r}}=N_{p\mathrm{r}}{\dot{z}}_{p\mathrm{r}}^{\mathrm{e}} \end{gathered}$$

where $k_{\mathrm{rb}}$ and $c_{\mathrm{rb}}$ are the equivalent stiffness and damping between the bridge rails, $z_{p\mathrm{r}}^{\mathrm{e}}$ is the displacement vector of the pth rail element corresponding to the wheel acting on the kth bridge element, and $z_{p\mathrm{r}}$ is the rail displacement value of the wheel action point.

Considering that there are m rail elements on a bridge element, and substituting Equation (8) into Equation (10), the following can be obtained:

$$M_{\mathrm{b}}^{\mathrm{e}}{\ddot{z}}_{\mathrm{b}}^{\mathrm{e}}+\left(C_{\mathrm{b}}^{\mathrm{e}}+C_{k\mathrm{b}}^{\mathrm{e}}\right){\dot{z}}_{\mathrm{b}}^{\mathrm{e}}+\left(K_{\mathrm{b}}^{\mathrm{e}}+K_{k\mathrm{b}}^{\mathrm{e}}\right)z_{\mathrm{b}}^{\mathrm{e}}-{\displaystyle \sum _{p=1}^m\left(C_{p\mathrm{rb}}^{\mathrm{e}}{\dot{z}}_{p\mathrm{r}}^{\mathrm{e}}+K_{p\mathrm{rb}}^{\mathrm{e}}{z}_{p\mathrm{r}}^{\mathrm{e}}\right)}=0$$

(9)

$$C_{p\mathrm{b}}^{\mathrm{e}}=c_{\mathrm{rb}}{{\displaystyle \int }}_{(p-1)l_r}^{p{l}_r}{N}_{\mathrm{b}}^{\mathrm{T}}{N}_{\mathrm{b}}\mathrm{d}x;p=1,2,\cdots m$$

$$C_{k\mathrm{b}}^{\mathrm{e}}={\displaystyle \sum _{p=1}^m{C}_{\mathrm{pb}}^{\mathrm{e}}}=c_{\mathrm{rb}}{{\displaystyle \int }}_0^{l_{\mathrm{b}}}{N}_{\mathrm{b}}^{\mathrm{T}}{N}_{\mathrm{b}}\mathrm{d}x$$

(10)

where $K_{k\mathrm{b}}^{\mathrm{e}}$ and $C_{k\mathrm{b}}^{\mathrm{e}}$ are the additional damping matrix and stiffness matrix of the kth bridge element due to elastic support, respectively. $C_{p\mathrm{rb}}^{\mathrm{e}}$ and $K_{p\mathrm{rb}}^{\mathrm{e}}$ are the coupling damping and stiffness matrix of the pth rail element of the kth bridge element due to elastic support, respectively. $l_{\mathrm{r}}$ and $l_{\mathrm{b}}$ are the length of the rail and bridge element, respectively.

The coupling damping matrix $C_{p\mathrm{rb}}^{\mathrm{e}}$ can be expressed as:

$$C_{p\mathrm{rb}}^{\mathrm{e}}=c_{\mathrm{rb}}{{\displaystyle \int }}_{(p-1)l_r}^{p{l}_r}{N}_{\mathrm{b}}^{\mathrm{T}}{N}_{\mathrm{pr}}\mathrm{d}x;p=1,2,\cdots m$$

(11)

2.2.3. Rail Element Equation

Similarly, the dynamic equation of the pth orbital element can be expressed as:

$$M_{p\mathrm{r}}^{\mathrm{e}}{\ddot{z}}_{p\mathrm{r}}^{\mathrm{e}}+\left(C_{p\mathrm{r}}^{\mathrm{e}}+C_{p\mathrm{rr}}^{\mathrm{e}}\right){\dot{z}}_{p\mathrm{r}}^{\mathrm{e}}+\left(K_{p\mathrm{r}}^{\mathrm{e}}+K_{p\mathrm{rr}}^{\mathrm{e}}\right)z_{\mathrm{r}}^{\mathrm{e}}-C_{p\mathrm{br}}^{\mathrm{e}}{\dot{z}}_{\mathrm{b}}^{\mathrm{e}}+K_{p\mathrm{br}}^{\mathrm{e}}{z}_{\mathrm{b}}^{\mathrm{e}}=-N_{x=x_{\mathrm{w}}}^{\mathrm{T}}F\left(t\right)$$

(12)

$$C_{prr}^{\mathrm{e}}=c_{\mathrm{rb}}{{\displaystyle \int }}_{(p-1)l_r}^{p{l}_r}{N}_{pr}^{\mathrm{T}}{N}_{pr}dx=c_{\mathrm{rb}}{{\displaystyle \int }}_0^{l_r}{N}_r^{\mathrm{T}}{N}_{r}dx$$

(13)

$$C_{p\mathrm{br}}^{\mathrm{e}}=c_{\mathrm{rb}}{{\displaystyle \int }}_{(p-1)l_r}^{p{l}_r}{N}_{p\mathrm{r}}^{\mathrm{T}}{N}_{\mathrm{b}}\mathrm{d}x;\hspace{0.33em}\hspace{0.33em}p=1,2,\cdots m$$

(14)

where $M_{p\mathrm{r}}^{\mathrm{e}}$, $C_{p\mathrm{r}}^{\mathrm{e}}$, and $K_{p\mathrm{r}}^{\mathrm{e}}$ are the mass, damping, and stiffness matrices of the pth rail element, respectively. $C_{prr}^{\mathrm{e}}$ and $K_{prr}^{\mathrm{e}}$ are the additional damping and stiffness matrix of the rail element, respectively, caused by the elastic support of the pth rail element in the kth bridge element.

Equation (15) can be obtained by substituting Equations (4)–(6) into Equation (12).

$$M_{p\mathrm{r}}^{\mathrm{e}}{\ddot{z}}_{p\mathrm{r}}^{\mathrm{e}}+\left(C_{p\mathrm{r}}^{\mathrm{e}}+C_{p\mathrm{rr}}^{\mathrm{e}}\right){\dot{z}}_{p\mathrm{r}}^{\mathrm{e}}+\left(K_{p\mathrm{r}}^{\mathrm{e}}+K_{p\mathrm{rr}}^{\mathrm{e}}\right)z_{\mathrm{r}}^{\mathrm{e}}-C_{p\mathrm{br}}^{\mathrm{e}}{\dot{z}}_{\mathrm{b}}^{\mathrm{e}}+K_{p\mathrm{br}}^{\mathrm{e}}{z}_{\mathrm{b}}^{\mathrm{e}}=-N_{x=x_{\mathrm{w}}}^{\mathrm{T}}F\left(t\right)$$

(15)

$$\begin{gathered} K_{p\mathrm{rw}}^{\mathrm{e}}=a_{\mathrm{in}}{k}_{\mathrm{H}}{N}_{x=x_{\mathrm{w}}}^{\mathrm{T}}{N}_{x=x_{\mathrm{w}}}^{} \\ {F}_{p\mathrm{rw}}^{\mathrm{e}}=a_{\mathrm{in}}{k}_{\mathrm{H}}r\left(x_{\mathrm{w}}\right)N_{x=x_{\mathrm{w}}}^{\mathrm{T}} \end{gathered}$$

where $K_{p\mathrm{rw}}^{\mathrm{e}}$ is the rail additional stiffness matrix caused by wheel–rail contact, and $F_{p\mathrm{rw}}^{\mathrm{e}}$ is the extra load vector caused by wheel–rail contact.

2.2.4. Wheel Dynamic Equation

The wheel is supported by contact force, and its dynamic equation can be expressed as:

$$m_{\mathrm{w}}{\ddot{y}}_{\mathrm{w}}+a_{\mathrm{in}}{k}_{\mathrm{H}}{y}_{\mathrm{w}}-a_{\mathrm{in}}{k}_{\mathrm{H}}{N}_{x=x_{\mathrm{w}}}{q}_{\mathrm{r}}^{\mathrm{e}}=a_{\mathrm{in}}{k}_{\mathrm{H}}r\left(x_{\mathrm{w}}\right)$$

(16)

where $m_{\mathrm{w}}$ is the quality of the wheelset.

2.2.5. Train–Rail–Bridge Coupling Time-Varying Equation

Accordingly, the coupled element dynamic Equation (17) of an arbitrary length combination, as shown in Figure 2, can be obtained by simultaneous solution of Equations (12)–(16).

$$\begin{array}{l}\left(\begin{array}{ccccc}{M}_{\mathrm{b}}^{\mathrm{e}}& & & & \mathrm{sym}.\\ 0& M_{1\mathrm{r}}^{\mathrm{e}}& & & \\ ⋮& ⋮& \ddots & & \\ 0& 0& \dots & M_{m\mathrm{r}}^{\mathrm{e}}& \\ 0& 0& \dots & 0& m_{\mathrm{w}}\end{array}\right)\left(\begin{array}{c}{\ddot{z}}_{\mathrm{b}}^{\mathrm{e}}\\ {\ddot{z}}_{1\mathrm{r}}^{\mathrm{e}}\\ ⋮\\ {\ddot{z}}_{m\mathrm{r}}^{\mathrm{e}}\\ {\ddot{y}}_{\mathrm{w}}\end{array}\right)+\left(\begin{array}{ccccc}{C}_{\mathrm{b}}^{\mathrm{e}}+C_{\mathrm{kb}}^{\mathrm{e}}& & & & \mathrm{sym}.\\ -C_{1\mathrm{br}}^{\mathrm{e}}& C_{1\mathrm{r}}^{\mathrm{e}}+C_{1\mathrm{rr}}^{\mathrm{e}}& & & \\ ⋮& ⋮& \ddots & & \\ -C_{m\mathrm{br}}^{\mathrm{e}}& 0& \dots & C_{m\mathrm{r}}^{\mathrm{e}}+C_{m\mathrm{rr}}^{\mathrm{e}}& \\ 0& 0& \cdots & 0& 0\end{array}\right)\left(\begin{array}{c}{\dot{z}}_{\mathrm{b}}^{\mathrm{e}}\\ {\dot{z}}_{1\mathrm{r}}^{\mathrm{e}}\\ ⋮\\ {\dot{z}}_{m\mathrm{r}}^{\mathrm{e}}\\ {\dot{y}}_{\mathrm{w}}\end{array}\right)+\\ \left(\begin{array}{ccccc}{K}_{\mathrm{b}}^{\mathrm{e}}+K_{\mathrm{kb}}^{\mathrm{e}}& & & & \mathrm{sym}.\\ -K_{1\mathrm{br}}^{\mathrm{e}}& K_{1\mathrm{r}}^{\mathrm{e}}+K_{1\mathrm{rr}}^{\mathrm{e}}+K_{1\mathrm{rw}}^{\mathrm{e}}& & & \\ ⋮& ⋮& \ddots & & \\ -K_{m\mathrm{br}}^{\mathrm{e}}& 0& \dots & K_{m\mathrm{r}}^{\mathrm{e}}+K_{m\mathrm{rr}}^{\mathrm{e}}+K_{m\mathrm{rw}}^{\mathrm{e}}& \\ 0& -a_{\mathrm{in}}{k}_{\mathrm{H}}{N}_{x=x_{\mathrm{w}}}& \dots & -a_{\mathrm{in}}{k}_{\mathrm{H}}{N}_{x=x_{\mathrm{w}}}& a_{\mathrm{in}}{k}_{\mathrm{H}}\end{array}\right)\left(\begin{array}{c}{\dot{z}}_{\mathrm{b}}^{\mathrm{e}}\\ {\dot{z}}_{1\mathrm{r}}^{\mathrm{e}}\\ ⋮\\ {\dot{z}}_{m\mathrm{r}}^{\mathrm{e}}\\ {\dot{y}}_{\mathrm{w}}\end{array}\right)=\left(\begin{array}{c}0\\ F_{1\mathrm{rw}}^{\mathrm{e}}\\ ⋮\\ F_{m\mathrm{rw}}^{\mathrm{e}}\\ a_{\mathrm{in}}{k}_{\mathrm{H}}r\left(x_{\mathrm{w}}\right)\end{array}\right)\end{array}$$

(17)

As previously stated, the train system is based on a theoretical model of a 10 DOFs multi-rigid-body system, which is illustrated in reference [39], along with its system equation, stiffness, and damping matrix.

According to the principle that the wheel–rail interaction unit and the total potential energy of the power system are fixed, the finite element motion equation of the train–rail–bridge interaction system is as follows [40]:

$$\left[\begin{array}{ccc}{M}_{\mathrm{bb}}& & \mathrm{sym}.\\ 0& M_{\mathrm{rr}}& \\ 0& 0& M_{\mathrm{vv}}\end{array}\right]\left\{\begin{array}{l}{\ddot{z}}_b\hfill \\ {\ddot{z}}_r\hfill \\ {\ddot{z}}_v\hfill \end{array}\right\}+\left[\begin{array}{ccc}{C}_{\mathrm{bb}}& & \mathrm{sym}.\\ C_{\mathrm{br}}& C_{\mathrm{rr}}& \\ 0& 0& C_{\mathrm{vv}}\end{array}\right]\left\{\begin{array}{l}{\dot{z}}_{\mathrm{b}}\hfill \\ {\dot{z}}_{\mathrm{r}}\hfill \\ {\dot{z}}_{\mathrm{v}}\hfill \end{array}\right\}+\left[\begin{array}{ccc}{K}_{\mathrm{bb}}& & \mathrm{sym}.\\ K_{\mathrm{br}}& K_{\mathrm{rr}}+K_{\mathrm{cr}}& \\ 0& K_{\mathrm{rv}}& K_{\mathrm{vv}}\end{array}\right]\left\{\begin{array}{l}{z}_{\mathrm{b}}\hfill \\ z_{\mathrm{r}}\hfill \\ z_{\mathrm{v}}\hfill \end{array}\right\}=\left\{\begin{array}{c}{F}_{\mathrm{b}}\\ F_{\mathrm{r}}+F_{\mathrm{cr}}\\ F_{\mathrm{v}}+F_{\mathrm{cv}}\end{array}\right\}$$

(18)

where M, C, and K represent mass, damping, and stiffness submatrix, respectively. $z$ and $F$ represent displacement and force subvector, respectively. In addition, subscripts b, r, v, and c represent the interaction effects caused by the bridge, rail, train, and wheel–rail contact.

2.2.6. Model Verification

According to the example provided by the literature [41], the solving model of the coupling element with an arbitrary length combination is verified. The acceleration time history response of the bridge is calculated using a group of 12 trains to pass through a six-span simply supported beam at a speed of 374 km/h. The bridge element is 3 m, the bridge damping ratio is 0.18, and the rail element is 0.6 m. Other bridge, rail, and train parameters can be found in the reference. The comparison between the calculated results of this paper and the results from the reference is shown in Figure 4, which shows that the amplitudes and trends of the two are basically the same, verifying the correctness of the method and program in this paper.

3. Damage Identification

Generally speaking, a bridge is invariably subjectively divided for structural damage identification, and the element’s damage is then preset [26,27,28]. On the one hand, if the length of the selected unit is too small, the search space will inevitably increase, and the corresponding computing time will also significantly increase. Moreover, due to the algorithm’s limitations, it may lead to the emergence of a locally optimal solution. On the other hand, if the length of the selected unit is too large, the recognition result is in an extensive range, which affects the accuracy of recognition. Moreover, it is generally assumed that the location of the damage happens to be an element in the FEM, and the damage location corresponds one-to-one to the element partition form of the FEM. However, this is an ideal and exceptional condition. In practice, there is an excellent possibility of a difference between the two, as shown in Figure 5.

3.1. Damage Model

The quality of the structure changes very little after structural damage, and the structural stiffness reduction factor $\mathsf{\theta }={[{\theta }_1,{\theta }_2,\cdots {\theta }_n]}^{\mathrm{T}}$ can reflect the degree of damage. The stiffness reduction factor is $\mathsf{\theta }={[0,0,\cdots 0]}^{\mathrm{T}}$ for complete structures. The following is the connection between the structural element stiffness matrix and the global stiffness matrix:

$$k_i^{\mathrm{b}}=\left(1-{\theta }_i\right)k_i^{\mathrm{u}}\left(0\le {\theta }_i\le 1\right)$$

(19)

$$K_i^{\mathrm{b}}=T_r^{\mathrm{T}}{k}_i^{\mathrm{b}}T$$

(20)

$$K_{\mathrm{b}}={\displaystyle \sum _{i=1}^n{\overline{K}}_i^{\mathrm{b}}}$$

(21)

where $k_i^{\mathrm{b}}$ and $k_i^{\mathrm{u}}$ stand for damaged and intact local stiffness matrices of the ith element in the structure. $T$ is the transformation matrix. $K_i^{\mathrm{b}}$ is an expanded ndof × ndof matrix denoted by ${\overline{K}}_i^{\mathrm{b}}$. ndof is the total number of structural freedom degrees. $K_{\mathrm{b}}$ represents the structural global stiffness matrix with damage, and n is the number of elements.

3.2. Fitness Function

The primary function of the identification factor is to reflect the location and degree of structural damage in a quantitative form. At present, damage can be effectively identified by using the correlation between two random variables as an identification factor, proven in reference [42,43].

• Cosine similarity.

  $$S_{\beta }={\displaystyle \sum _{i=1}^N\frac{{\mathit{\Psi }}_i^{\mathrm{T}}{\mathit{\Phi }}_i}{\left({\mathit{\Psi }}_i^{\mathrm{T}}{\mathit{\Psi }}_i\right)\left({\mathit{\Phi }}_i^{\mathrm{T}}{\mathit{\Phi }}_i\right)}}$$

  (22)

  where ${\mathit{\Psi }}_i$ and ${\mathit{\Phi }}_i$ are the measured and theoretical data of the acceleration at the ith measuring point, respectively, with a total of N measuring points.
• Pearson correlation coefficient.

  $$S_{\rho }={\displaystyle \sum _{i=1}^N\frac{Cov\left({\mathit{\Psi }}_i,{\mathit{\Phi }}_i\right)}{\sqrt{D\left({\mathit{\Psi }}_i\right)}\sqrt{D\left({\mathit{\Phi }}_i\right)}}}$$

  (23)

  where $Cov$ represents covariance, and $D$ represents variance.
• Fitness function.

  $$OBJ=\alpha S_{\beta }+\left(1-\alpha \right)S_{\rho }$$

  (24)

  where $\alpha$ is the weight coefficient, take 0.4~0.6.

3.3. Flow Chart

In order to solve the above problems, a larger unit length is used to divide the bridge based on the dichotomy idea. Then, a more considerable damage step and damage range are selected for rough identification to lock the scope of the bridge damage location. On this basis, a length of 1/2 is used to subdivide these possible damage locations, the appropriate damage parameters are selected, and the GA is used to identify and repeat the process until the desired accuracy is achieved. In this way, the increase in search space and the exponential increase in computing costs caused by a direct subdivision of the unit can be avoided. Combined with the improved GA, the calculation process is shown in Figure 6.

4. Illustrative Example

Taking the 30 m simply supported beam as an example, assuming that the irregularity of the rail collected in the field is shown in Figure 7, eight trains pass through the simply supported beam at a constant speed v, and the schematic diagram of the whole system is shown in Figure 8. N_S sensors with a sampling frequency of 2000 Hz are uniformly arranged on the bridge. The IGA was used for 500 independent repeated calculations and analyses. The parameters of the bridges, rails, and trains used in the calculation are acquired from the literature [41]. The parameters of the genetic algorithm are shown in

Table 1. Parameters of the IGA.

Parameter	Value	Parameter	Value
Coding mode	Integer encoding	Maximum iteration	70
Crossover rate	0.8	Catastrophe threshold	25
Mutation rate	0.1	Number of catastrophes	1
Initial population size	50	Identification factor weight coefficient α	0.6

.

The measured data are simulated by adding a few Gaussian white noises to the calculated value of the FEM in a particular damage state [44], such as Equation (25):

$${\mathit{\Psi }}_i={\mathit{\Phi }}_i+e_p{N}_{o}D\left({\mathit{\Phi }}_i\right)$$

(25)

where $e_p$ is the noise degree, $N_o$ is a standard normal distribution vector with zero mean and unit standard deviation, and D represents variance. The signal-to-noise ratio (SNR) of the acceleration signal of the jth measuring point is calculated according to Equation (26):

$$SN{R}_j=10\times \lg \left(\frac{1}{m}{\displaystyle \sum _{i=1}^m{y}_{ji}^2}/\frac{1}{m}{\displaystyle \sum _{i=1}^m{\sigma }_{ji}^2}\right)$$

(26)

where y_ji is the acceleration value of the ith time point of the jth measuring point, σ_ji is the noise value of the ith time point of the jth measuring point, and m is the total number of time points.

4.1. Recognition Process

4.1.1. First Recognition

As mentioned in Section 3, if the length of the selected unit is too small, the search space will inevitably increase, and the corresponding computing time will also significantly increase. This effect will be discussed in Section 4.1.2. Therefore, the FEM is divided into 20 elements, each 1.5 m long. Figure 9 shows a comparison between the location of the actual damage and the unit number of the FEM (A’s position is just staggered from the FEM, and the length is 1.5 m. B’s position is less than half of the element, and the length is 0.6 m.). The specific damage parameters are shown in

Table 2. Damage information parameters under conditions 1~2.

Condition	Step Size and Range/%	Damage Location	Damage Degree/%	SNR/dB	NS	v/(km/h)
1	[0:5:10]	A	8	14	5	40
2		B	10	14	5	40

.

ANS indicates the number of iterations required to find the optimal solution for the first time in each calculation. Figure 10 shows the frequency histogram of ANS. As shown in Figure 10, most ANS are concentrated between 8 and 14. Due to the Preserving the Optimal Strategy, the algorithm converges quickly. The iterative process of one of the conditions 1~2 is analyzed, as shown in Figure 11.

As shown in Figure 11, around the 20th generation, the population has undergone a catastrophe. Better individuals are produced, the diversity of the population is improved, and the optimal local solution of the algorithm is avoided.

Figure 12 shows the overall distribution of the first recognition result of conditions 1~2. It can be concluded that the actual damage location is invariably near the recognition result or its symmetrical location. As can be seen from Figure 12a, 14.4% of the results are 5% damage to unit a, and 5% damage to unit s. It can be recorded as “14.4% 『a#5∪s#5』”, the same below. The distribution of the nearby and the symmetrical position of the a and s units is shown in Figure 13.

As shown in Figure 13, in the identification result 『a#5∪s#5』, a unit itself contains the left half of the real damage location A, and the symmetrical position b of the s element contains the right half of the real damage location A. The rest of the case is the same. The above conclusions can be attributed to the following reasons:

• The boundary conditions of simply supported beams are symmetrical.
• The FEM division of the bridge is not consistent with the actual damage.
• The randomness of Gaussian white noise.

4.1.2. Second Recognition

Only one of the recognition results (14.4% 『a#5∪s#5』and 4.4% 『i#5』) in conditions 1~2 is analyzed because of the above conclusion (the actual damage location is invariably near the recognition result or its symmetrical location), and the analysis process of other cases is similar. Figure 14 shows the comparison of the damage location and the element division. As described in Section 3, if the length of the selected unit is too large, the recognition result is in an extensive range, which affects the accuracy of recognition. Therefore, the bridge elements are re-divided, and the range of variables is limited. Take the first identification result 14.4% 『a#5∪s#5』of condition 1 as an example. The first identification result shows unit a damaged, and the variable is limited to the range of [a b s t] in the second recognition, corresponding to the second subdivision unit, which is [1 2 3 4 23 24 25 26]. Similarly, the first identification result shows unit s damaged, and the variable is limited to the range of [a b c r s t] in the second recognition, corresponding to the second subdivision unit, which is [1 2 3 4 5 6 21 22 23 24 25 26]. The total recognition rate of identifying damage locations twice is obtained according to Equation (27).

$$I{R}_{\mathrm{Z}}=I{R}_{\mathrm{F}}\times I{R}_{\mathrm{S}}$$

(27)

where IR_Z is the total recognition rate, IR_F is the probability that the real damage is included in the range of the nearby and the symmetrical position of the first recognition result, and IR_S is the probability of the second recognition.

From Figure 15, the overall locating effect of the two conditions is better, and the recognition rate is greater than 95%. However, in the quantization results, there is a slight error between the recognition results in Figure 15b and the actual damage, and most of them are identified as unit 13 (length 0.75 m) with 8% uniform damage. The actual damage is contained in unit 13, with a length of 0.6 m and damage of 10%. However, this error is relatively small and is acceptable for practical application.

In terms of computational efficiency, because [2:2:10] is set as the recognition step and range of the second recognition, the search space for condition 1 is 2400. Plus, the search space in the first recognition is 3600, a total of 6000. If the traditional one-time recognition method is used, and the first recognition is divided according to a smaller length, the search space is as high as 40,000, which is 6.7 times that of this method. Thus, the search space of this method is much lower than that of the traditional one-time recognition method. Moreover, with the increase in the number of units damaged, the number of search spaces in the traditional one-time identification method exponentially increases; this method shows significant superiority, and the computing time is significantly reduced.

4.2. Noise Impact Analysis

The method is based on the FEM, but the FEM inevitably has errors. In order to analyze whether the algorithm is sensitive to noise, only the SNR value in conditions 1~2 is changed to 12 dB, and the other conditions remain unchanged. The two recognition results are as follows.

As shown in Figure 16, when the SNR decreases, a few first recognition results do not meet the above conclusions (The part marked with red circles). In addition, the IR_Z in Figure 17 also decreases, indicating that noise impacts the recognition rate. After changing the SNR value in conditions 1~2, its recognition effect can be calculated and fit.

As shown in Figure 18, there is an “S” curve relationship between SNR and IR_Z. The higher the SNR, the higher the IR_Z. When the SNR is greater than 14 dB, the recognition effect is better and has good robustness. The smaller the SNR, the stronger the interference of the noise, and the recognition effect will gradually decrease.

4.3. Analysis of the Influence of the Number of Measuring Points

The more measuring points, the more data will be obtained. However, the more measuring points, the more sensors need to be arranged, resulting in a waste of resources. In order to give the suggested value of a reasonable number of measuring points, by the same token, one can change the number of measuring points in conditions 1~2 in Table 2 and study the recognition effect of the algorithm.

As shown in Figure 19, there is a positive correlation between N_S and IR_Z, and a better recognition effect can be achieved by setting five measuring points in the whole bridge. When it is fewer than five, it will drop rapidly. This phenomenon happens because the fewer measuring points there are, the fewer data sets are used for analysis, and the interference of noise is relatively strong, thus reducing the recognition rate of damage.

4.4. Analysis of the Influence of Train Speed

The sampling frequency of the sensor is fixed. If the speed is too fast and the number of samples is reduced, it may lead to a significant error in the recognition results. One can change the v in conditions 1~2 in Table 2, study the recognition effect of the algorithm, and give a reasonable suggested value.

As shown in Figure 20, v is negatively correlated with IR_Z. When v > 75 km/h, the recognition effect will gradually decrease with speed increase. This phenomenon happens because for fixed-span bridges, the faster the moving speeds of the test train, the less data can be collected at each measuring point, and the interference of noise is relatively substantial, thus reducing the identification rate of damage.

5. Conclusions

In order to study the influence on the recognition effect when there is a difference between the element division form of the FEM and the actual damage location, a multi-times location damage method based on variable-length elements and an improved genetic algorithm is proposed. First of all, the selection, crossover, and mutation operators in the traditional genetic algorithm are improved, and catastrophe and neighborhood search mechanisms are introduced to enhance the local optimization ability of the local algorithm. The test function verifies the performance of the improved genetic algorithm. Then, the time-varying equation of train–rail–bridge coupling is established, and the objective function is constructed using the bridge’s dynamic response under trainload. Finally, an example is used to perform multiple independent repeated calculations for different damage conditions, and the stability and convergence rate of the algorithm is statistically analyzed. Moreover, the effects of noise, the number of measuring points, and train speed on the identification results are discussed. From this study, the following conclusions can be drawn:

(1)

The sum of the search space of this algorithm is less than that of the traditional one-time recognition method, and the computational efficiency is improved. With the increase in the number of damaged units, the number of the former slowly increases, while the number of the latter exponentially increases. As far as an example is concerned, the traditional method is 6.7 times that of the algorithm.

(2)

The algorithm converges quickly. The number of iterations required for the first appearance of the optimal solution ANS is mostly between 8 and 14, which is on the low side. Under the action of the optimal holding strategy operator, the algorithm can converge quickly.

(3)

There is a specific relationship between noise, number of measuring points, train speed, and recognition rate. The recognition rate has an “S” curve relationship with noise, a positive correlation with the number of measuring points, and a negative correlation with train speed.

(4)

Even if there is a difference between the damage location and the element division form of the FEM, the damage location can be correctly located, but the quantitative result of the damage degree will be affected. Under the given suggested value, the recognition effect of the algorithm is good and can reach more than 95%.