This section presents the proposed framework of transfer function fitting. Figure 1 shows a conceptual diagram of the proposed framework, which minimizes the error between the target transfer function obtained from observation records and the estimated transfer function obtained from the analytical model. To reproduce the target transfer function efficiently, a process is employed in which errors between the target and estimation are minimized only at the dominant response. Furthermore, by extracting multiple modes from a single observation point and normalized errors between the target and estimated response, vibration characteristics up to higher-order modes can be considered even with a limited number of observation points. Note that normalization is necessary to increase the weight for higher order modes. The minimized error vector {e ₀} is defined in Eqs 1a–1c: e 0 = R f 1 R a 1 R f 2 R a 2 … R f i R a i … R f N R a N T , (1a) R f i = R f i 1 … R f i j … R f i L , R a i = R a i 1 … R a i j … R a i L , (1b) R f i j = f i j T − f i j E f i j T , R a i j = a i j T − a i j E a i j E , (1c) where f i j T and a i j T are the dominant frequencies and amplitudes of the ith mode in the target transfer function, respectively. The transfer function was calculated by utilizing the Fourier spectral ratio of the absolute acceleration on the jth floor of the observed building to the input acceleration. f i j E and a i j E indicate the dominant frequency and amplitude, respectively, for the estimated transfer function of the ith mode for the jth floor obtained by the seismic response analysis. Therefore, {e ₀} comprises the relative errors of the dominant frequencies and amplitudes between the target and estimated values normalized by the target values for each mode. The dominant modes were visually selected with respect to the target transfer function. Moreover, as shown in Figure 1, by setting a trapezoidal band-pass filter defined by f ₁–f ₄ that envelopes each mode, the target and estimated transfer functions compare the dominant modes in the same bands.

This section presents an overview of the parameter identification flow using the MIEC method, according to the framework of a previous study (Uesaka et al., 2021b). For detailed steps of inverse analysis using MIEC, refer to previous research (e.g., Suzuki, 2018; 2019). The calculation in the MIEC method considers a nonlinear transform system F that outputs an N-dimensional vector {a ₀} from an input M-dimensional parameter vector {p} based on Eq. 2: a 0 N × 1 = F p M × 1 , (2) where all vectors are treated as column vectors, F is assumed to be a frequency response analysis using a linear MDOF system SR model in our study. {p} comprises the shear stiffness and damping ratio of the superstructure and the spring stiffness, damping ratio, and damping coefficient of the dynamic soil spring (details are described later). {e ₀} in Eq. 1a is used as the output vector {a ₀} in this study. In other words, {p} is searched for such that the error vector is 0.

To correct an error in {a ₀} iteratively, the perturbation impulse matrix B , comprising incremental output vectors when the perturbation {Δp} is added to the input, is recursively calculated. Subsequently, vector {Δα}, which satisfies Eq. 3, is calculated: r N × 1 = a target N × 1 − a 0 N × 1 = B ′ N × M Δ α M × 1 , (3) where B ′ is obtained by selecting the significant singular vectors in B and the low-rank approximation of the matrix. Δ α is the incremental ratio vector of the perturbation that corrects the residual error of output vector {r}, calculated using the low-rank approximated general-inverse matrix B ′ + , as shown Eq. 4: Δ α M × 1 = B ′ + M × N r N × 1 , (4)

Subsequently, {p} is iteratively modified using {Δα}, as shown Eq. 5, until the residual sum of squares (RSS) error (the norm of {r}) becomes smaller than an allowable tolerance: p new M × 1 = p M × 1 + ∆ p M × 1 ∘ ∆ α M × 1 , (5) where “◦” represents the Hadamard product, {Δp}◦{Δα} is the correction vector of {p}, and {p _new} denotes the modified input vector.

This section presents the target and estimated transfer function of the equations of motion of the SR model used in this study. The equation of motion of the SR model is expressed using Eqs 6–9. This equation is generally expressed as a relative displacement system with respect to the ground input motion displacement y, such as the free-soil surface response, as shown in Figure 2A, where k _i and c _i denote the shear stiffness and damping coefficients of the superstructure, respectively. K _H , C _H , K _R , C _R , K _HR(RH), and C _HR(RH) are the stiffness and damping coefficients of the sway,- rocking,-, and SR-coupled springs, respectively. m _i is the mass of each mass point, J ₀ is the rotational inertial mass concentrated at the foundation, and H _i is the story height. Although the influence of SR-coupled springs has been proven to be non-negligible in pile foundations (Cristina et al., 2013), the effect is generally ignored. We show that parameter identification can be performed considering the effect of an SR-coupled spring. Note that the rotation angle of the seismic isolation layer and rocking angle of the foundation due to the rocking soil spring may be different, which is not shown in Figure 2. To simplify the following formulations, we adopted the equation of motion of the SR model, which ignores these differences in rotation angles. In the numerical experiments presented in Section 3, where the rotation angle of the seismic isolation layer and rocking angle of the foundation are different, we confirmed that the physical parameters can be identified according to the formulation in Section 2. The equation of motion for the SR model is as follows, where the time function is abbreviated as x(t) = x. In the numerical analysis in this research, the equation of motion was handled in the frequency domain; however, following equations are expressed in the time domain, because the actual observation records are time history data: M x ¨ + C x ˙ + K x = − M e y ¨ , (6) M = d i a g m i ⋯ m 1 m 0 J 0 , (7) C = c i − c i − c i ⋱ ⋯ − c i H i ⋮ ⋮ − c i H i ⋯ c 1 + c 2 − c 1 c 2 H 2 − c 1 H 1 − c 1 c 1 + C H c 1 H 1 + C H R c 2 H 2 − c 1 H 1 c 1 H 1 + C R H C R + ∑ j = 1 n c j H j 2 , (8) K = k i − k i − k i ⋱ ⋯ − k i H i ⋮ ⋮ − k i H i ⋯ k 1 + k 2 − k 1 k 2 H 2 − k 1 H 1 − k 1 k 1 + K H k 1 H 1 + K H R k 2 H 2 − k 1 H 1 k 1 H 1 + K R H K R + ∑ j = 1 n k j H j 2 , (9) where [M], [C], and [K] represent the mass, damping, and stiffness matrices, respectively. { x = x i x i − 1 ⋯ x 0 θ 0 T is the relative displacement vector from the ground motion input point, and e = 1 ⋯ 1 0 T is a state vector representing the inertial force distribution. In Eq. 6, the influence of foundation embedding is assumed to be relatively small, and the ground acceleration in the rocking direction is assumed to be zero.

In the previous studies, ground responses were often required to estimate building and soil spring parameters (e.g., Cristina et al., 2013). However, ground responses are not always observed. Therefore, in this study, the equations of motion were transformed to estimate these parameters without using the ground response. According to the method presented by Tojo and Nakai (2022) and Tojo et al. (2023), the equation of motion of the relative displacement system was considered with respect to the plane rotated by θ ₀ from the horizontal plane shown in Figure 2B. Additionally, to reflect the influence of the SR coupling spring, we attempted to move the forces associated with the K _HR and C _HR terms on the right side, as shown in Eqs 10, 11; that is, they can be considered separately in the two equations. One includes the damping and stiffness matrices or coefficients [C _S ], [K _S ], C _R , and K _R for the superstructure and rocking soil spring. The other contains coefficients C _H and K _H of the horizontal soil spring. In addition, Eq. 11 becomes Eq. 12 when considering a vibration state x ¨ 0 ≫ y ¨ in which the excitation force, such as the external wind force on the superstructure, is more dominant than the input ground motion: M ¯ X ¨ + C ¯ X ˙ + K ¯ X = − M ¯ e X ¨ 0 + H θ ¨ 0 + f ¯ x 0 , (10) m 0 X ¨ 0 + C H x ˙ 0 + K H x 0 = c 1 X ˙ 1 + k 1 X 1 + f ¯ θ 0 , (11) m 0 x ¨ 0 + C H x ˙ 0 + K H x 0 = c 1 X ˙ 1 + k 1 X 1 + f ¯ θ 0 , (12) M ¯ = M S M B S J 0 , (13a) M S = d i a g m i ⋯ m 1 , (13b) M B S = m i ∑ H i ⋯ m 1 H 1 , (13c) C ¯ = C S C R , C S = c i − c i − c i ⋱ c 1 + c 2 , (14) K ¯ = K S K R , K S = k i − k i − k i ⋱ k 1 + k 2 , (15) f ¯ x 0 = − 0 C H R x ˙ 0 − 0 K H R x 0 , (16a) f ¯ θ 0 = − C H R θ ˙ 0 − K H R θ 0 , (16b) where M ¯ , C ¯ , and K ¯ are the matrices of Eqs 13a–13c, Eq. 14, and Eq. 15, X = X i X i − 1 ⋯ X 1 θ 0 T is the relative displacement vector corresponding to the building shear deformation excluding the rocking displacement from the foundation mass point, and H = ∑ H i ∑ H i − 1 ⋯ H 1 T is the height vector.

Subsequently, considering the transfer function of the response to the input when the external forces f ¯ x 0 and f ¯ θ 0 defined by Eqs 16a, b on the right sides of Eqs 10, 12 are ignored, Eqs 17a, 17b, and Eq. 18 are obtained, respectively. In this study, Eqs 17a, 17b were used as the target or estimated transfer function to identify the physical parameters related to the superstructure and rocking soil spring. Eq. 18 was utilized to identify the parameters of the horizontal soil spring. The effect of the SR coupling spring could be considered because f ¯ x 0 and f ¯ θ 0 were treated as arbitrary external forces that satisfied the target transfer function of Eqs 17a, 17b, and Eq. 18: G k FL = X ¨ k + X ¨ 0 + ∑ H k θ ¨ 0 X ¨ 0 + ∑ H k θ ¨ 0 = U ¨ k X ¨ 0 + ∑ H k θ ¨ 0 , (17a) G R = θ ¨ 0 X ¨ 0 + H e q , R θ ¨ 0 , (17b) G H = x ¨ 0 X 1 , (18) H e q , R = H 1 2 + … + ∑ H i 2 H 1 + … + ∑ H i , (19) where G _kFL and G _R represent the transfer functions when the horizontal absolute acceleration response of kth floor of the superstructure or the rocking absolute acceleration response of the foundation mass point is divided by the horizontal effective acceleration input. The horizontal effective acceleration input is defined as the input obtained by multiplying the foundation rocking input times the height of each floor and adding it to the ground motion. G _H is the transfer function of the horizontal absolute acceleration response of the foundation mass point to the horizontal relative displacement response of first floor to the foundation mass point. ΣH _k is the height up to kth floor, and H _eq,R defined by Eq. 19 is the equivalent height of the rigid-body rocking mode proposed by Tojo and Nakai (2022). By considering the transfer functions of two separate equations in this way, the superstructure and dynamic soil spring can be estimated without using the ground response. However, Eqs 17a, 17b must handle the response when an earthquake is input, and Eq. 18 must handle the response when an external force is applied to the superstructure.

The analysis object was an SR model of a five-mass system with mass i = 5 in Figure 2. Tables 1, 2 list the specifications of the superstructure and soil springs of the SR model. Numerical analyses were performed in the frequency domain so that the analytical model for the superstructure and SR soil springs was linear. The specifications were set as a model of the short-side direction based on the stiffness and mass in the seismic design of an existing reinforce concrete (RC) BI building with three stories. Section 4 describes the building details. Because the study subject was a shear-type building with a seismic-resistant wall above the BI layer, it was modeled as an equivalent shear-type beam element with a sufficiently large bending stiffness. The rotational inertial masses of the floors above the BI layer were consolidated into first floor (referred to 1FL in Table 1). The damping constant was the average value based on the measured database of RC buildings (Architectural Institute of Japan, 2020). The BI layer was modeled as a bending shear beam, considering the shear and bending stiffnesses and damping constant based on the design values. The damping of the superstructure is treated as the complex stiffness k _i ^* = k _i (1 + i2h _i ). The foundation comprised cast-in-place concrete piles. Tables 3, 4 list the foundation and soil specifications, respectively. The damping constants of each soil layer were set based on the relationship between the Q-value and the S-wave velocity (V _S ) (Architectural Institute of Japan, 1987). Each soil spring of Table 2 was calculated by using the thin-layer element method (TLM) (e.g., Tajimi and Shimomura, 1976) based on the specifications in Tables 3, 4. Furthermore, the spring stiffness and damping of the soil springs were approximated as constants according to the method of the Architectural Institute of Japan (2006). In addition, in Table 2, the values of the SR springs assuming that the piles were damaged are listed in each lower row. These SR springs considering damaged pile were calculated by setting the cross-sectional area of more than half of the pile heads with length 0.25 m to 1/10. The arrangement of the piles is the same position as the seismic isolation device shown in section 4.

Figures 3A, B show the acceleration time–history waveforms and acceleration response spectra of the input seismic motions. Three earthquake motions were used as input seismic motions: the first was a white noise wave with little bias due to its spectral characteristics (hereinafter referred to as WN); the second was an east to west direction component of the KiK-net Haga observation wave published by National Research Institute for Earth Science and Disaster Resilience, 2019 obtained from the 2011 Tohoku Earthquake (hereinafter referred to as Haga); and the third was artificial designed seismic motion provided by the Building Center of Japan, 1992 (hereinafter referred to as Dgn). The duration and sampling time were 81.92 and 0.01 s, respectively. In the numerical analysis, scenarios receiving the seismic input or external force of wind were assumed and the following cases were considered. Note that, in the discussion presented in Section 3.4, the transfer functions of Eqs 17a, 17b were evaluated in Analysis (1), and that of Eq. 18 in Analysis (2). Specifically, Analysis (1) was a seismic response analysis, where each seismic motion was input through the soil spring in the horizontal direction, and Analysis (2) was a excitation analysis of the superstructure, where the acceleration excitation of the WN was input in the horizontal direction to the mass point of 1FL of the superstructure.

Recent studies (e.g., Architectural Institute of Japan, 2020; Cruz and Miranda, 2021) have shown different trends in the damping properties of architectural structures. In one trend, the attenuation of higher-order modes increased with increasing frequency. Another trend maintained approximately constant attenuation of the lower and higher modes. Considering these findings, to examine the applicability of the proposed framework to indirect monitoring of soil and foundations, we conducted an analysis in which the stiffness was varied assuming damage to the piles. Furthermore, we confirmed the effects of differences in the damping characteristics of the ground on the identification results.

Table 5 lists the numerical analysis cases. Four cases were evaluated in this study, denoted by the initials “N” and “D” for cases with undamaged and damaged piles, respectively, as listed in Table 2. In addition, “-C” denotes the case in which the damping constants (h _H , h _R , h _HR ) of the dynamic soil spring in the horizontal and rocking directions were used to provide a constant divergent damping to the vibration frequency. “-P” indicates the case in which the damping was proportional to the vibration frequency using each damping coefficient (C _H , C _R , C _HR ). Table 5 shows the natural frequencies and damping constants obtained through the complex eigenvalue analysis (e.g., Terazawa and Takeuchi, 2018) based on the complex stiffness of the SR model in each case, along with the results of the foundation-fixing condition (FB). Focusing on the differences between the damping models, although the natural frequencies of N-C and N-P or D-C and D-P were nearly the same, differences in the damping characteristics of the primary and secondary modes were observed.

The input parameters identified, denoted as {p} in Eq. 2, were the stiffness and damping constants or damping coefficients (k _i, K _H , K _R , K _HR , h _i , h _H , h _R , and h _HR or C _H , C _R , and C _HR ) of each part shown in Table 1. We assumed that the bending stiffness k _B of the BI layer was known. In the MIEC method, the initial values and perturbation quantities of each parameter to be identified should be provided; Table 6 lists the initial values and perturbation quantities of each parameter. The initial values provided in the inverse analysis were rounded to a power of 10, assuming that the approximate value of each parameter could be estimated in advance. For convenience, the amount of perturbation of the shear stiffness was defined relative to the shear stiffness k _i multiplied times the story height H _i . {p _new}, obtained from Eq. 5, was set as the difference from the initial value. This difference was added to the initial values and input into the analytical model. Because each parameter (K _HR , h _HR , and C _HR ) related to the SR-coupled spring was an external force satisfying the relation in Eqs 17a, 17b, and Eq. 18, negative values were allowed. The other parameters were treated as absolute values. The perturbation {Δp} for each parameter was approximately 1/10 to 1 times the initial value of each parameter. Assuming an observation error for each response acceleration time history obtained in the analysis, the target transfer function was calculated by adding noise with a signal-to-noise (SN) ratio of 100. The noise was generated based on a normal distribution with a mean of 0 and standard deviation of 0.1. The SN ratio was defined as the ratio of the variance of each response in the numerical analysis to the variance of the noise. The number of modes to be selected for matrix B ′ in Eq. 3 was determined using the norm ratio of the error vector in the study by Suzuki (2019), and the threshold of the norm ratio was 0.3. The convergence threshold of the RSS error was 0.03, and the maximum number of repetitions was 150. However, if a tendency for convergence to a smaller error was observed, the maximum number of repetitions was increased to 300. The maximum dominant response extracted by the band-pass filter was set to be up to the fifth-order mode, which could be determined visually. The response of the observation point used for identification was assumed to be BFL, 1FL, and RFL only, considering the installation position of the sensor in the actual building, as shown in Section 4. The target transfer function was calculated according to Eqs 17a, 17b, and Eq. 18 using the respective horizontal and rotational acceleration responses or the displacement response obtained from those responses. Figures 4A–C respectively show the target transfer function and trapezoidal filter for the N-P, N-C, and D-P cases as examples of a band-pass filter setup for extracting the dominant response. The estimated transfer functions are also shown in these figures, as are the transfer functions for the unmeasured floors. Each physical parameter was identified using the following procedure.

Step 1

All physical parameters were temporarily identified based on the transfer function of Eqs 17a, 17b. (The target and estimated values are shown in Figures 4A–C as G _RFL, G _1FL, and G _R , respectively). In this step, repeated calculations were performed based on the initial values and amount of perturbation in Table 6. The dominant responses from six points were extracted: the first- and second-order peaks for G _RFL; the first-, third-, and fourth-order peaks for G _1FL; and the fifth-order peak for G _R in Figures 4A–C, respectively. Basically, the number of points for extracting the dominant response of the transfer function could be the same as the maximum objective mode order. However, in the first mode of a BI building, the superstructure was considered to move as an almost rigid body. To incorporate this constraint, we selected the dominant mode for both 1FL and RFL when extracting the first-order mode (e.g., see the notation “first” in Figure 4A).

Step 2

K _H , K _HR , h _H h _HR , and h _R or C _HR , C _HR , and C _R were identified again based on Eq. 18. The transfer functions based on Eq. 18 are shown in Figures 4A–C as G _H . The initial values and perturbations of the parameters to be identified in Step 2 were used from the values in Table 6. The initial values of the other parameters were fixed using the values obtained in Step 1. The dominant modes selected were up to two peaks. One was the mode in which the horizontal response of foundation was predominant (represented as Sway in Figure 4A). The other was excited by the rocking response due to the effect of the SR coupling spring (represented as Rocking in Figures 4A–C). However, for N-P and D-P, which had large dissipative damping at high frequencies, the high-frequency peak at approximately 20 Hz present in the N-C case was not observed and amplitude increased monotonically. Therefore, by adding amplitudes at 7.5, 20, and 30 Hz as fitting point of the target transfer function, we attempted to improve the parameter identification accuracy by capturing the shape of the transfer functions (see the symbols ▾ in Figures 4A–C).

Figures 7A, B show plan views of BFL, 1FL, and RFL and the elevation in the direction of the short side of the building, respectively. The locations of the accelerometers installed on each floor are indicated in the plans. Measurements were not taken for 2FL or 3FL. This building is a three-story RC BI building completed in 1987. The foundation is a pile foundation comprising cast-in-place concrete piles. The BI layer comprises 14 natural-rubber-based laminated rubbers and eight viscous dampers. The arrangement of the piles is the same position as the laminated rubbers. Approximately 10, 22, and 30 years after the completion of construction, one laminated rubber piece was removed, and an aging investigation was conducted (Higashino et al., 1997; Hamaguchi et al., 2009; Kamoshita et al., 2018). According to Kamoshita et al. (2018), the vertical stiffness of laminated rubber increased with aging, with a fluctuation range of approximately 19%–2%. Regarding the horizontal stiffness of the same member, the effect depended on the evaluation section of shear strain γ. It tended to decrease from the time of shipment in the γ = 10%–50% section, but increased from the time of shipment in the γ = 10%–200% section. The former decreased by −16%–7%, and the latter increased by 7%–14%. The aging of viscous dampers for approximately 20–30 years has also been investigated (Hamaguchi and Sasaki, 2009; Tanaka et al., 2018). The shear resistance obtained by sampling a part of the viscous body was estimated to be approximately −2%–4% of the standard value. These investigation results confirmed that no remarkable aging occurred, although minor variations were observed. Figure 7C shows the number of seismic waves, maximum absolute acceleration, inter-story velocity, and inter-story deformation angle recorded annually from June 2007 to March 2023 in BFL when the physical parameters were evaluated. For the inter-story velocity and inter-story deformation angle, the relative velocity and relative displacement were obtained by integrating the absolute acceleration in BFL and 1FL once or twice, respectively, and the inter-story deformation angle was calculated by dividing the relative displacement by the story height of the BFL. The magnitude of the responses increases significantly in 2011, as recorded during the 2011 Tohoku Earthquake. Excluding these responses, the maximum responses recorded from 2007 to 2023 vary slightly; however, most are small earthquakes with maximum accelerations of less than 100 cm/s². The maximum inter-story deformation angle is less than 1% rad excluding that in 2011, and the deformation of the BI layer is extremely small.

The physical parameters for the observation records were identified based on the five-mass point system SR model used in Section 3. The observational records evaluated were on the short side of the building. The identified physical parameters were the same as those in Section 3. Table 6 was also referenced for the amount of perturbation of each physical parameter. The initial values of the parameters based on the design values in Tables 1, 2 were used. However, each initial value of the SR-coupled spring was set to zero in accordance with normal seismic design. As obtaining the damping characteristics of soil springs in advance is difficult, the damping model of the dynamic soil spring was examined in advance using two types of models: one with constant dissipative damping (damping constant) with respect to frequency (case C in Section 3) and another with damping proportional to the frequency (case P). As a result, the model with case C was used owing to its low RSS error. These analysis conditions were common for all earthquake motions.

Wind observations were not conducted at the target buildings, and only strong motion observation records were obtained. Therefore, only Eqs 17a, 17b were considered in this analysis. The target transfer functions were G _R , G_1FL, and G_RFL from the responses of BFL, 1FL, and RFL, where the accelerometers were installed. Among the physical parameters, the identification values of h _R , K _H , and h _H in addition to the stiffnesses K _HR and h _HR of the SR-coupled spring were not analyzed in detail. Regarding the acceleration response used for the input and output, the horizontal direction was the average value of the horizontal accelerometers at both ends in the X-direction, as shown in Figure 7A. The rocking direction responses were calculated from the difference in the vertical accelerometers at both ends in the Y-direction in the BFL divided by the distance. The number of dominant response modes extracted was of the fifth order. The band-pass filter settings for extracting the dominant response were set by obtaining the average value of the transfer function based on the seismic observation records and visually judging the dominant modes. Figures 8A, B show the target transfer function based on the observation record and the estimated transfer function based on the identification results, respectively. The estimated transfer function was calculated using the regression equation of the physical parameters and the relationship between the average amplitude, which will be described later. The extraction target was from the first to fifth mode. The primary mode shown in Figure 8A was applied to G _1FL and G _RFL to consider the constraint of the rigid-body mode, as described in Section 3. The fourth and fifth modes had large response variations and were unclear. Therefore, the extraction range of the filter was relatively large. Determining the extraction method for the dominant response in such cases requires further investigation. Comparing the transfer functions in Figures 8A, B, the dominant responses of the target and estimated transfer functions generally correspond well for the first-to third-order modes (see approximately 1, 7, and 11 Hz in the figures). There is room for improvement in the consistency of higher-order modes beyond the fourth order and the accuracy of the rocking response of G _R .

This section presents an analysis of the aging and amplitude dependence of the physical parameters identified based on the analysis conditions described in the preceding paragraph.

Figure 9A shows the average value μ, maximum value (Max.), and minimum value (Min.) for all annual seismic waves of the maximum inter-story deformation angle of the foundation floor (BFL). However, for 2011, the results for March 11, including the 2011 Tohoku Earthquake, and the succeeding results were calculated separately. The 2010 and 2012 periods, before and after 2011, are hatched in yellow and shown separately. The results in Figures 9B–G are also hatched in the same way. The average maximum inter-story deformation angle shown in Figure 9A is approximately 0.1% rad or less in each year, except on 11 March 2011. No significant differences are observed; however, after 2015, earthquakes with relatively small responses dominate.

The left side of Figure 9B shows the RSS error of the target and estimated transfer function in each year. The legend shows the mean μ and standard deviation ±σ within a single year. The right side of Figure 9B shows the relationship between the average maximum inter-story deformation angle R _max shown in Figure 9A and the average RSS error. The solid black lines in Figure 9B represent the lines when the errors are 20%, 30%, and 35%. The average error is approximately 30% in each year. The bias of the identification result depending on the spectral characteristics of the ground motions was smoothed to a certain degree through averaging. Observing the amplitude dependence diagram on the right, these errors are approximately the same, regardless of the amplitude.

Figures 9C–G show the identification results of the shear stiffness (k _i ) of superstructures 3F, 2FL, 1FL, and BFL and the rocking soil spring stiffness K _R in each year on the left side. The right sides of Figures 9C–G show the relationship between the identification results and the average maximum inter-story deformation angle shown in Figure 9A. From the numerical test described in Section 3 (see Table 9), the smaller the RSS error, the better the identification accuracy of physical parameters. Therefore, the weighted average value μ ∼ and standard deviation σ ∼ defined in Eq. 20a, 20b and 21a, 21b were adopted for the identification results of each parameter shown in Figures 9C–G. The solid black line in each figure represents the initial value used for identification. The amplitude dependencies in Figures 9C–G can be expressed by curves based on regression analysis. Figure 9F, where the amplitude dependence is relatively clear, utilizes a power function approximation. The remaining figures use log-linear approximations. μ ∼ = ∑ l = 1 N w l p l ∑ l = 1 N w l , (20a) w l = 1 / e l R S S , (20b) σ ∼ = ∑ l = 1 N w l δ l 2 N − 1 ∑ l = 1 N w l , (21a) δ l = p l − μ ∼ , (21b) where p _l is the physical parameter identified with the ground motion l in a given year, and w _l is the weighting factor for the ground motion, which is the inverse of the RSS error e l R S S at the time of identification.

The secular change and amplitude dependence of each physical parameter shown in Figures 9C–G were confirmed. The largest earthquake was the 2011 Tohoku Earthquake, and damage to this building due to this earthquake was not confirmed by the previous report. Therefore, we will first focus on the identification results for 2010–2012. In Figures 9C–G, no decrease in stiffness is observed in any part. From these results, it can be inferred that no significant damage was caused to any story during the 2011 Tohoku Earthquake. The decrease in the stiffness k _BFL of the BI layer in Figure 9F in 2011 can be attributed to the amplitude dependence, which will be discussed later. The stiffness increases slightly in the other stories; however, these increases are likely due to instability in the identification accuracy, such as variation of the spectral characteristics of the input seismic motion. These instabilities are thought to affect the fluctuations in other years as well. Therefore, the secular change in stiffness of each layer other than the seismic isolation layer was judged not to be a significant result. Next, we will discuss the amplitude dependence of the seismic isolation layer in Figure 9F. The relationship in Figure 9F demonstrates that the stiffness decreases as the inter-story deformation angle increases. Therefore, changes in the seismic isolation layer over time are not due to structural deterioration, but rather are caused by the stiffness identification results increasing or decreasing depending on the average amplitude of the earthquake ground motion in that year. Morii et al. (2015) reported the strain dependence of high-damping laminated rubber during small deformations, and the same effect was assumed for the laminated rubber of the building in this study. In Figure 9G, the rocking spring stiffness K _R is approximately 20×10⁸ kNm/rad. However, K _R increases and decreases significantly every year. Moreover, no consistent trend in amplitude changes can be confirmed. This result can be said to indicate the accuracy limitations of the proposed method of extracting a single response of the transfer function, especially because the dominant response in the transfer function of G _R in Figure 8A is unclear.

Figures 10A–D show the variation trends in the damping of the superstructure with respect to aging and vibration amplitude in the same manner as in the respective images in Figure 9, starting from the upper layer. Similar to Figure 9, Figures 10A–C depict log-linear approximations, and Figure 10D, in which the amplitude dependence is relatively clear, shows the curve and regression equation expressed using power function approximation. Figure 10E presents the regression results obtained using a log-linear approximation for each physical parameter excluded from the discussion of identification results. To address the identification results in Figure 10E, other data such as wind observations must be obtained, which is a topic for future work.

Figures 10A–C reveal no significant changes in the damping characteristics of superstructure attributed to aging. Furthermore, no apparent variation in the damping constant due to amplitude dependence is observed. In years with large fluctuations, the variance of the identified parameters is also large, suggesting a decline in estimation accuracy. In addition, the variation remains at approximately 1%–2%, consistent with the statistical value from the existing RC structure database (Architectural Institute of Japan, 2020). The damping constant h _BFL of the BI layer shown in Figure 10D decreases as the average maximum inter-story deformation angle increases, in the same manner of its stiffness. This amplitude dependence is considered to be the cause of yearly fluctuations in the attenuation constant. The mechanism through which the attenuation constant decreases depending on the amplitude is considered to be as follows. The standard value of a damping force F _D of the viscous damper of the BI building is proportional to the 0.4–1.0 power of the velocity v, which is expressed in Eqs 22a, b (Tanaka et al., 2018). However, because the analytical model used for identification is a linear model in the frequency domain, the damping force F ∼ D is expressed by Eq. 23. Subsequently, by contrasting Eqs 22a, b), and (Eq. 23) with the relationship between the damping constant and damping coefficient, we obtain Eq. 24. According to this equation, the damping constant decreases as the amplitude increases. This velocity amplitude-dependent characteristic can be presumed also to have affected the damping constant shown in Figure 10D. F D = α e − 0.043 t 0.01 v / 10 β = C v β , (22a) β = 1.0   v < 0.1   cm / s 0.59   0.1  cm / s ≤ v < 1.0   cm / s 0.4   1.0  cm / s ≤ v , (22b) F ∼ D = C ∼ v , (23) h = C ∼ / C c = C v β − 1 / C c , (24) where α is a constant depending on the velocity v, t is the temperature, β is a power exponent depending on v, and C _c is the critical damping coefficient.

The above results demonstrate that no significant decrease in stiffness or change in attenuation occurs due to aging in each physical parameter, excluding the BI layer, and no damage due to earthquake motion is estimated. Moreover, no clear tendency is observed in the amplitude dependence of these parameters. Considering the results of the numerical experiments in Section 3, the fluctuations other than those of the seismic isolation layer are thought to be due to the instability of the identification accuracy associated with the spectral characteristics of seismic motion. For the BI layer, clear amplitude dependence in the stiffness and damping constant was confirmed. The tendency for the stiffness to decrease with amplitude corresponds well with the previous research results, such as those for high-damping rubber. On the contrary, the tendency for the damping to decrease depending on the amplitude is the opposite of the amplitude dependence caused by high-damping rubber. This finding is inferred to be due to the properties of the viscous damper installed in the building. As the amplitude dependence of seismically isolated buildings differs depending on the seismic isolation system adopted, as much data as possible should be collected in the future.

No clear aging or amplitude dependence was observed for the rocking soil spring either. However, the ambiguity of the dominant response of the target transfer function likely affected the identification results, especially for the rocking soil spring. This aspect represents a limitation of the accuracy of the single-point fitting of the transfer function proposed in this study. Furthermore, the process used to extract the dominant reactions visually in the proposed framework should be improved.

Finally, a seismic response analysis of representative seismic waves from 2007 to 2023 using the physical parameters identified using the mass point system model was performed, and the effectiveness of the parameter evaluation results was confirmed by comparing the responses with the observation records. Considering the results presented in the previous section, the physical parameters were evaluated using the regression equations of the amplitude dependence shown in Figures 9, 10. The amplitude adopted to the regression equation to estimate each physical parameter was determined from the maximum inter-story deformation angle of each earthquake motion. In the seismic response analysis, the horizontal acceleration response in the observation record of the BFL mass point was used as the forced acceleration. The earthquake motions to be conducted were the following six waves, which had different amplitude levels: the 2011 Tohoku Earthquake and the aftershocks on the same day (the maximum inter-story deformation angles of BFL, R _max = 5.56% and 1.93% rad, were defined as the mainshock and aftershock, respectively); the 2008 and 2022 earthquakes, in which the inter-story deformation angle of the BFL was maximum between 2007 and 2010 or 2012 and 2023 (R _max = 0.27% and 0.88% rad for BFL, respectively); and the 2007 and 2023 earthquakes with small amplitude levels (R _max = 0.05% and 0.02% rad for BFL, respectively), totaling six waves.

Figures 11A–F show the maximum response distribution along the height direction based on the analysis performed using the identified physical parameters. The upper row depicts the absolute acceleration, and the lower presents the displacement relative to the BFL. The results based on the initial values (Ini.) and observation records (Obs.) are also provided in the same images. Figures 11A–F demonstrate that the maximum acceleration and relative displacement are consistent with the observation records. The consistency of the responses is also higher than that of the initial values. The initial stiffness of the BI layer based on the seismic design was assumed to be 40% of the shear strain, and the strain level of the evaluated ground motions was lower than that assumed in the design. Therefore, the identification results considering the amplitude dependence exhibited better consistency.

Figures 12A–F present the acceleration response spectra (h = 5%) in the horizontal directions of RFL and 1FL and the rocking direction of BFL obtained from observation (Obs.) and analysis (Identification parameter: Ide. or initial value: Ini.) for six representative waves, arranged chronologically from 2007 onward. The responses of RFL and 1FL in Figures 12A–F show that the dominant responses around 1.0–2.0 s agree well between the analysis based on identified parameters and observation results. In Figures 12A, E, F, the response around 0.1–0.2 s, which is the second-order mode, also agrees between the observation and analysis results. The rocking responses of BFL in Figures 12A–F demonstrate that higher-order peaks more than 0.1 s were reproduced in the analysis based on the identification results, whereas the amplitude was underestimated.

The above results confirm that the response reproducibility was improved by identifying physical parameters considering higher-order modes using the proposed method. In addition, the regression results for physical parameters based on the assumed amplitude dependence described in the previous section can explain building responses to earthquake motions of different amplitude levels.