Application of Story-Wise Shear Building Identification Method to Actual Ambient Vibration

A sophisticated and smart story stiffness System Identification (SI) method for a shear building model is applied to a full-scale building frame subjected to micro-tremors. The advantageous and novel feature is that not only the modal parameters, such as natural frequencies and damping ratios, but also the physical model parameters, such as story stiffnesses and damping coefficients, can be identified using micro-tremors. While the building responses to earthquake ground motions are necessary in the previous SI method, it is shown in this paper that the micro-tremor measurements in a full-scale 5 story building frame can be used for identification within the same framework. The SI using micro-tremor measurements leads to the enhanced usability of the previously proposed story-wise shear building identification method. The degree of ARX models and cut-off frequencies of band-pass filter are determined to derive reliable results.


Introduction
Much research interest has been directed recently to system identification (SI) of civil, mechanical, and aerospace structures in response to the increasing need of enhancement of safety and upgrade of ability of damage detection (or damage diagnosis) of various kinds of structures (Hart and Yao, 1977;Beck and Jennings, 1980;Hoshiya and Saito, 1984;Kozin and Natke, 1986;Agbabian et al., 1991;Koh et al., 1991;Ghanem and Shinozuka, 1995;Hjelmstad et al., 1995;Shinozuka and Ghanem, 1995;Doebling et al., 1996;Hjelmstad, 1996;Masri et al., 1996;Housner et al., 1997;Kobori et al., 1998;Johnson and Smyth, 2006;Nagarajaiah and Basu, 2009;Fujino et al., 2010;Ji et al., 2011). Such need results from the accelerated demand of rapid assessment on material aging issues and of continuing use of buildings after earthquakes. The Building Continuity Plan (BCP) and resilience-oriented design framework also supports the increasing need of advancement of structural health monitoring technologies. It is also well-recognized that SI plays an important and principal role in reducing gaps between the constructed structural systems and their structural design models (model refinement).
Basically, two kinds of branches of the SI technique are known. One is called the modal parameter SI, which is well-established, and another is called the physical parameter SI. The modal parameter SI is well-established and various types of research have been conducted and accumulated so far [for example, see Hart and Yao (1977) and Beck and Jennings (1980)]. Since the modal parameters are global parameters and stable in principle, many useful techniques have been proposed. In the modal parameter SI, for evaluating the natural frequency and damping ratio, the observations at two places at least (usually the base and the top) are necessary. Furthermore, for evaluating modal shape identification, the simultaneous observations (or interpolation from fewer observations) at many places are usually required that is often a cumbersome task in the modal parameter SI method. When the observations at all the floors of a building are possible, all story stiffnesses of the shear building model with known floor masses may be identified within an acceptable accuracy (Hjelmstad et al., 1995).
In contrast to such modal parameter SI, the physical parameter SI has been developed for direct identification of physical parameters (stiffness and damping coefficients). For example, Takewaki and Nakamura (2000;2005) introduced a method based on the work by Udwadia et al. (1978). In that method, a shear building model is used and stiffness and damping coefficients of a specific story are identified directly from the floor accelerations just above and below the specific story. However, the method by Takewaki and Nakamura (2000;2005) has a difficulty resulting from the small signal/noise (SN) ratio in the low frequency range and cannot be applied to ambient vibration data, e.g., microtremors, and to high-rise buildings with large height-to-width aspect ratios. The former problem has been a major and most difficult problem in the field of SI and the latter problem has been tackled by Minami et al. (2013) and Fujita et al. (2013) by using a shear-bending model. Hernandez-Garcia et al. (2010a; have developed an interesting method of damage detection using a floor-by-floor approach to enhance the efficiency and accuracy of the identification results. Xing and Mita (2012) have devised a time-domain substructure damage identification method for shear buildings by focusing on a substructure consisting of one story. Zhang and Johnson (2013) have developed another substructure identification method for shear buildings, which considers the noise effect of recorded data on identification accuracy and utilizes an iterative inductive procedure from the top story. Furthermore, a combined method of the modal parameter SI and the physical parameter SI is also well used (for example, Shinozuka and Ghanem, 1995;Barroso and Rodriguez, 2004). After modal parameters are identified, physical parameters are determined by solving inverse problems in which the existence and uniqueness of solutions are principal problems. Takewaki and Nakamura (2009) have developed a method for identifying the temporal variation of modal properties of a base-isolated building during an earthquake by using the relation between auto-regressive eXogenous (ARX) model parameters and modal parameters. The difficulty arising in the limit manipulation in the method by Takewaki and Nakamura (2000;2005) has been overcome by introducing an ARX model in the previous paper (Maeda et al., 2011). The weakness of a small signal-to-noise (SN) ratio in the low frequency range in the method (Takewaki and Nakamura, 2000;2005) has been avoided by using the ARX model and introducing new constraints on the ARX parameters. Another difficulty due to small vibration levels of micro-tremor has been tackled by introducing a combination of the ARX model, filtering in the frequency domain (low and high-cut filter), and averaging in the time domain (sequential time-window shift for Fourier transformation and averaging). In this paper, for obtaining reliable and stable identification of story stiffnesses, a practical procedure to decide the combination of the degree of ARX models and cut-off frequencies of filtering is investigated by applying a multiobjective optimization algorithm. It is shown that the previously proposed story-wise stiffness identification method is applicable in the new approaches to actual recorded data of an ambient vibration level with the help of this sophisticated combination.

Governing Equations of Shear Building Model
Consider an N-story shear building model with viscous damping as shown in Figure 1. Let m j and k j denote the mass of the j-th floor and the story stiffness of the j-th story and let c j be the viscous damping coefficient of the j-th story. Since the formulation in the frequency domain is appropriate in the present formulation, all the governing equations are expressed in the frequency domain. Let "i" denote the imaginary unit.
Keeping the relationsU(ω) = iωU(ω), Ü(ω) = −ω 2 U(ω) in mind, the equations of motion in the frequency domain for this shear building model subjected to the horizontal ground acceleration üg is expressed as: where U(ω) and Üg(ω) are the Fourier transforms of the horizontal displacements u(t) = {u j (t)} of floors relative to ground and the ground base acceleration üg, respectively. The vector 1 in the righthand side of Eq. (1) denotes 1 = {1, 1, . . ., 1} T . The mass, stiffness, and damping matrices of shear building model are defined by: The transfer function T(ω) can be defined from Eq. (1) as the ratio of the horizontal displacement to the ground acceleration.

Sub-system Identification Based on Identification Function
Based on the mathematical formulation for the SI of the physical parameters (Takewaki and Nakamura, 2000), the story stiffness and damping coefficient can be expressed by using the limit manipulation of the identification function f j (ω) for ω → 0. The story stiffness and damping coefficient are derived as: where M j = ∑ N i=j m i . Referring Takewaki and Nakamura (2000), the identification function can be defined by: In Eqs (4-6), Ü 0 = 0. From Eqs (4) and (5), it is concluded that the story stiffness and damping coefficient can be derived only from the accelerations just above and below floors of the object story in the previously proposed identification method (Takewaki and Nakamura, 2000).

Identification of Stiffness and Damping Coefficient using ARX Model Taylor Series Expansion of Transfer Function using ARX Parameters
In the previous SI method (Takewaki and Nakamura, 2000;2005) for shear building models as formulated in Eqs (4) and (5), the limit manipulation of the identification function f j (ω) for ω → 0 was needed. However, when the identification functions are evaluated from the raw data such as actual vibration testing data, it is often the case that the identification functions become unstable and exhibits a large variability in the low frequency range. To overcome this difficulty, an ARX model is introduced, which is a time-domain model. The reliability of the ARX model in this direction has been confirmed and the applicability of the ARX model to shear building models and shear-bending models has been demonstrated in Maeda et al. (2011) and Minami et al. (2013).
The identification function f j (ω) defined by Eq. (6) can be rewritten by introducing the transfer function G j,j−1 (ω) between j-th and (j − 1)-th floors as where G j,j−1 (ω) ≡ (Üg + Ü j )/(Üg + Ü j−1 ). The reason of the unstable phenomenon of the identification function in the low frequency range can be understood from Eq. (7). From the theoretical investigation, the limit value of the transfer function G j,j−1 (ω) at ω = 0 should be G j,j−1 (0) = 1, because the j-th floor and (j − 1)-th floor move identically at ω → 0. Therefore, the limit value of the denominator of the identification function 1/G j,j − 1 (ω) − 1 for ω → 0 is zero. Furthermore, the limit value of the numerator of the identification function ω 2 M j for ω → 0 is also zero. For these reasons, the limit value of the identification function f j (ω) for ω → 0 is theoretically indefinite. By using ARX parameters, the transfer function can also be expressed as (see Appendix) For evaluation of the limit value of the transfer function G j,j − 1 (ω) for ω → 0, the formulation of the Taylor series expansion of the transfer function G j,j−1 (ω) is meaningful. The Taylor series expansion of the transfer function G j,j − 1 (ω) in terms of ARX parameters can be defined as Considering the relationship of ARX parameters in Eq. (8) and the coefficients of the Taylor series expansion in Eq. (9), the coefficients A 0 , A 1 , A 2 of the Taylor series expansion can be formulated in terms of the ARX parameters {a k }, {b k } as follows.
where asum ≡ ∑ n k=1 a k , bsum ≡ ∑ n k=1 b k . By introducing the real and imaginary parts of A j as A j = A R j + iA I j and substituting the properties of the real and imaginary parts of A j into Eqs (10-12), the transfer function can be reduced to:

Constraints on ARX Parameters Derived from the Limit Value of Taylor Series Expansion of Transfer Function
As mentioned in the previous section, it is meaningful to note that Eqs (14a,b) can be derived from the mechanical interpretation, i.e., the j-th floor and ( j − 1)-th floor move identically at ω → 0 and G j,j − 1 (ω) should not include linear terms of ω judging from Eqs (4-6).
For enhancing the reliability of the proposed SI method, these relations will be used as the constraints in the estimation of the ARX parameters. It is found that Eq. (15) is a linear equation for the ARX parameters. A batch processing least-squares estimation method (see Appendix) (Adachi, 2009) provides: R, θ, and f are defined in Appendix. By applying the Lagrange multiplier method, the linear constraint can be incorporated into the batch processing least-squares estimation method as p T θ = −1 where p = {1, . . ., 1, −1, . . ., −1} T . Therefore, the present method is reduced to the problem for solving the following equations.
[ R p

Objective Functions for Determination of Order of ARX Model
It is important to investigate how to determine the order of the ARX model. Let us introduce the following objective function in the time domain as the ensemble average of amplitude vector in time domain.
where ft(t i ) is defined by using an estimated time history u ARX from ARX parameters and raw measurement data uraw as: In Eq. (20), σraw is the SD of uraw, and uraw is the mean value of the time history. The first term [u ARX (t i ) − uraw(t i )]/σraw of the objective function in time domain indicates the accuracy of u ARX compared with uraw. The second term [uraw(t i ) −ūraw]/σraw can be regarded as the weighting function for the amplitude level of the time history, i.e., the influence of data on the objective function is small when uraw(t i ) ≈ūraw ≈ 0. When the estimated time history u ARX , which can be computed from the ARX parameters as u ARX = θ T φ(k), matches with the measured data in an acceptable accuracy, the objective function Eq. (19) is close to 0.
On the other hand, let us introduce another objective function in the frequency domain.
where Graw denotes the transfer function computed from the records after appropriate post-data processing, e.g., filtering and ensemble averaging, and G ARX is the transfer function evaluated by the ARX model. Note that ω 1 = ω l and ω NN = ωu, where ω l and

Ambient Vibration Measurement in Five-Story Steel Frame Structure
Actual micro-tremor observations were conducted in an experimental real-size building at the Uji campus, Disaster Prevention Research Institute of Kyoto University, Japan. An overview photograph of the building is shown in Figure 2. The servo-type velocity meters (VSE-15D; Tokyo Sokushin) were installed in several stories. The velocity resolution of this velocity meter is 10 −4 mm/s. Figure 3 shows a photo of an installed velocity meter. The frame dimension and its shear building model are shown in Figure 4. The locations of velocity meters for three patterns of identification are presented in Table 1. In all measurement patterns, the velocity meters are fixed at the basement and roof floor so as to obtain the transfer function of the objective frame structure. Pattern A is aimed at identifying the story stiffnesses of the first and second stories and pattern B the fourth and fifth stories. On the other hand, pattern C is set for identifying the third story. The ambient measurements data were recorded in the long-span direction and short-span direction, respectively. The floor masses are estimated as m 1 = 28.7 × 10 3 kg, m 2 = 28.2 × 10 3 kg, m 3 = m 4 = 27.6 × 10 3 kg, and m 5 = 26.6 × 10 3 kg by preliminary investigation. Figure 5 shows some examples of recorded ambient data at the base, 3rd floor, 4th floor, and roof in the short-span direction in Pattern C. Durations of the measurement time were 5 and 10 min. Figure 6 presents the transfer function of the roof velocity to the base for each pattern of the short-span direction measurements. In Figure 6, the transfer functions evaluated by raw data, i.e., no post-processing, are compared with those after smoothing and filtering. From these preliminary analyses of steel structure, natural frequencies can be estimated as shown in Table 2.

Determination of Identification Parameters by Iterative Manipulation
As seen in Figure 6, the transfer functions evaluated by raw data exhibit large variability in the low frequency range. For this reason, it is difficult to obtain the limit value of the identification function directly from the raw data. It is well known that the order of ARX models is important to obtain reliable identification results. Furthermore, it is also observed that the band-pass filter as postdata processing stabilizes the variability of the ARX modeling. Taking these unknown parameters for the structural identification into account, it is needed to investigate the influence of the combination of these unknown parameters on the stability of the identification results. One of the solutions to this problem, i.e., the identification method combining the previous method (Hernandez-Garcia et al., 2010b;Xing and Mita, 2012) with the band-pass filtering and the ARX modeling is presented here. Figure 7 shows the outline of the identification procedure including the band-pass filtering and the ARX modeling. In the flowchart (Figure 7), three procedures are needed; Procedure 1 is for determination of the order of the ARX model, which was derived by the parametrical analysis to minimize objective functions defined by Eqs (19) and (21); procedure 2 is for determination of the lower cut-off frequency; and procedure 3 is for that of the upper cut-off frequency.
Finally, the story stiffnesses are evaluated from the limit value of the identification function in Eq. (4) by selecting unknown identification parameters. If the stability of the limit value is not acceptable, it is needed to select the order of ARX models again in procedure 1. Figure 8 presents examples of SI by using the measured ambient data of the real structure. The red lines, which can be derived by iterative manipulations, are the finally selected identification parameters for the order of the ARX model and cut-off frequencies. However, as seen in Figure 8, the objective functions in terms of these identification parameters vary drastically. Therefore, these procedures of iterative flows to determine several identification parameters may depend on the knowledge of structural engineers so as to derive reliable identification results.

Determination of Identification Parameters by Multi-Objective Optimization
It is important to select an appropriate combination of identification parameters for reliable identification of the story stiffnesses. In this section, this combination will be derived by applying a multi-objective optimization solver based on genetic algorithm. The design parameters for optimization are the order of the ARX model, lower and upper bounds of cut-off frequency. The objective functions are given by Eqs (18) and (20). In the application to the measured ambient data as shown in the following section, four velocity sensors are used simultaneously, which means that the transfer functions and identification functions can be evaluated  at three different consecutive stories. For the application to this experimental data, the number of the objective functions in the multi-objective optimization algorithm is six, i.e., the error rate in time domain and frequency domain at three consecutive stories.
The multi-objective optimization algorithm (MOGA-II) implemented in modeFRONTIER 4.5 (ESTECO) is used to investigate the optimal combination of the identification parameters. Initial design parameters are provided by the uniform random number method. Figure 9 describes the multi-dimensional analysis in the post-data processing for all individuals in the multi-objective optimization (the number of designs is 8000). The points connected by lines represent a generated design and the corresponding error rates of objective functions. From this figure, an excellent combination of the identification parameters can be derived by determining an allowable range of objective functions. Figure 9 shows two results (top and bottom illustrations) for different allowable ranges of the objective functions. In the bottom illustration, the error rate in time domain is widen, which causes the wide variability of design parameters. As for the order of ARX models, two different ranges of ARX order can be seen as 26-36 and 58-60 by widening the error rate in time domain. From this multi-dimensional analysis, the evaluation of the objective functions influences the decision making on the identification parameters. Figure 10 shows the comparison of the transfer function (roof/base) obtained from raw data with the transfer functions derived by the optimal combination of the order of ARX models and the cut-off frequencies in the multiobjective optimization (nine combinations are selected here). In Figure 10, no difference can be seen in these combinations, which demonstrates the reliability of the evaluation of error rate in the multi-objective optimization procedure.

Stability of Stiffness Identification by Limit Value of Identification Function
By applying the multi-objective optimization algorithm, several designs of the combination of identification parameters are obtained, which can make the objective functions in time and frequency-domain minimized in an acceptable manner. In this section, stability of the story stiffnesses is investigated by applying the SI method to the various (and optimized) identification functions derived by the multi-objective optimization. As explained before, the present SI method needs the evaluation of the limit value of identification functions. In the previously proposed SI method in the same framework, this limit value was just determined by selecting the raw value of the identification function at ω ≈ 0 (called Method 1). In this sense, Method 1 is a simple and primitive one. Figure 11 shows the cumulative frequency of the identified story stiffness k i (i = 1, 2, . . ., 5) for various combinations derived from the multi-dimensional analysis in Figure 9 where the number of optimal designs is about 150. It can be observed that identified story stiffnesses by Method 1 have large variability although the selected combination groups of the identification functions are similar. This means instability of the story stiffness identification.
For enhancing the stability of story stiffness identification, other approaches to limit value evaluation are proposed as: Method 2: Using gradient sensitivity of identification function.
Method 3: Using model pole derived by ARX model.
A conceptual illustration of these methods for limit value evaluation is shown in Figure 12. In Method 2, the gradient sensitivity of the identification function is evaluated successively in the frequency domain. The limit of the identification function is determined by evaluating the mean value of the identification functions in a stable range where the gradient of the identification function is sufficiently small (see Figure12: Method 2). While, in Method 3, a raw value of the identification function at a particular frequency determined by the frequency for the model pole of ARX model is used as the limit value of the identification function. The particular frequency point is given by assuming an appropriate constant value, e.g., 0.5, as the ratio to the model pole frequency (see Figure12: Method 3). Let us see Figure 11 again where the cumulative frequency of story stiffnesses is evaluated by the proposed methods. The variability of the identified story stiffnesses can be evaluated from the comparison of the gradient in the cumulative frequency shown in Figure 11. As for the first, second, and third stories, identified story stiffnesses have smaller variability compared with those by Method 1. However, improvement of stability for the fourth and fifth story stiffnesses cannot be observed even in applying Method 2 and 3. This means that it is difficult to determine story stiffnesses in these cases. This instability needs further investigation on the relationship between the objective functions, Eqs (19) and (21), and the identification  function derived by the ARX model. Tables 3 and 4 summarize the identified story stiffnesses and natural frequencies obtained by the shear model using the identified story stiffnesses. The identified story stiffnesses derived by Method 2, 3 in Table 3 are determined by the mean values of the selected optimal designs in Figure 11. In Table 3, a frame analysis was conducted by using a frame model including measured section size of columns and beams. It should be mentioned that the boundary condition of the numerical frame analysis may not correspond to the real structure. It is difficult to define "reference stiffness values" because a real building structure cannot be expressed by a shear building model exactly.

Conclusion
A previously proposed story-wise stiffness identification method using an ARX model for a shear building structure has been applied to the case where the shear building is subjected to micro-tremors measured in a full-scale steel frame structure. While earthquake ground motions and the building responses to such inputs are necessary in the previous method, it has been shown that micro-tremors can be used for identification within the same framework. This enhanced the usability of the previously proposed identification method. The difficulty in the selection of the combination of identification parameters, i.e., the order of the ARX model and the cutoff frequencies of band-pass filtering, has been tackled by a multi-objective optimization algorithm. Advantageous features are as follows: (1) Micro-tremors can be used as an input for the previously proposed SI method for shear building models regardless of its small vibration level.
(2) In order to investigate identification parameters for reliable identification based on the ambient vibration, the objective functions in time and frequency domains have been proposed by applying filtering in the frequency range (low and high cuts) and averaging in time domain (sequential time-window shift for Fourier transformation and averaging).
(3) Stability in the story stiffness evaluation for optimal combinations of identification parameters derived by the multiobjective optimization has been compared among various limit-value estimation approaches. It has been observed that the limit value by selecting the raw value of the identification function at ω ≈ 0 is not necessarily appropriate. Proposed approaches such as using gradient sensitivity of the identification function or the model pole of the ARX model make it stable to derive story stiffnesses.