# Mathematical Modeling for Lateral Displacement Induced by Wind Velocity Using Monitoring Data Obtained from Main Girder of Sutong Cable-Stayed Bridge

Authors: Gao-Xin Wang; You-Liang Ding
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Based on the health monitoring system installed on the main span of Sutong Cable-Stayed Bridge, GPS displacement and wind field are real-time monitored and analyzed. According to analytical results, apparent nonlinear correlation with certain discreteness exists between lateral static girder displacement and lateral static wind velocity; thus time series of lateral static girder displacement are decomposed into nonlinear correlation term and discreteness term, nonlinear correlation term of which is mathematically modeled by third-order Fourier series with intervention of lateral static wind velocity and discreteness term of which is mathematically modeled by the combined models of ARMA(
7,4
) and EGARCH(
2,1
). Additionally, stable power spectrum density exists in time series of lateral dynamic girder displacement, which can be well described by the fourth-order Gaussian series; thus time series of lateral dynamic girder displacement are mathematically modeled by harmonic superposition function. By comparison and verification between simulative and monitoring lateral girder displacements from September 1 to September 3, the presented mathematical models are effective to simulate time series of lateral girder displacement from main girder of Sutong Cable-Stayed Bridge.


---

## Body

## 1. Introduction

Nowadays, long span cable-stayed and suspension bridge structures are commonly constructed at home and abroad. On account of their flexible structural characteristics, displacement response from main girder of long-span bridge structure swings obviously impacted by strong aerostatic and fluctuating wind actions. According to aerostatic response analysis on Sutong Cable-Stayed Bridge by Xu et al., the lateral displacement response from main girder can approach 1.2 m under strong wind velocity 40 m/s with attack angle 0° [1]; and research results from buffeting response analysis on Golden Gate Bridge by Vincent showed that extreme buffeting amplitude from main girder can reach 1.7 m under strong wind velocity 31 m/s [2]. Such large amplitude can definitely threaten comfort and safety of the whole bridge structure. For example, severe wind vibration from main girder of Tacoma Suspension Bridge in Washington state eventually brought about collapse of the whole bridge structure under wind velocity 19 m/s [3]. Therefore, it is of great significance to research displacement response impacted by wind loads from main girder of long span bridge structures, and especially the lateral displacement response, as one fairly important part for main girder, should be specifically valued.Theoretical exploration, numerical simulation, and wind tunnel tests for lateral displacement response have been carried out to some extent. Cheng and Xiao improved the calculation method for aerostatic stability and further concluded that instable lateral displacement was 4.24 m under critical static wind [4]; Long et al. analyzed lateral displacement response from Sidu Suspension Bridge through ANSYS finite element simulation and concluded that maximum lateral displacement at middle span was 32.26 cm, which is close to 1/1000 length of main span [5]; Yu et al. researched lateral displacement response from Xihoumen Suspension Bridge through wind tunnel tests and concluded that lateral displacement at horizontal angle 10° was larger than that of other angles [6].However, for mechanism complexity of lateral displacement response impacted by aerostatic and fluctuating wind actions, traditional methods of theoretical deduction, numerical simulation, and wind tunnel tests are difficult to accurately reflect the actual lateral displacement response of bridge structure, on account of uncertain boundary condition, imprecise assignment of initial parameters, and inappropriate ignorance of subordinate factors. In recent years, with development of structural health monitoring technology, it is feasible to install monitoring sensors on long span bridge structures, monitoring data of which can authentically reflect bridge structural behaviors under actual environment and load actions. Although wind field of long span bridge structures has been widely monitored in recent years [7–9], lateral displacement response is rarely monitored and researched; thus real correlation regularity between lateral displacement response and wind action is still covered. Additionally, lateral displacement response under actual operation environment is also affected by other random factors, which have never been taken into account by researchers before. Therefore, lateral displacement response from main girder is necessarily researched upon monitoring data to reveal real structural behavior of long span bridges.In this paper, based on the health monitoring system installed on the main span of Sutong Cable-Stayed Bridge, GPS displacement and wind field are real-time monitored and analyzed. According to analytical results, apparent nonlinear correlation with certain discreteness exists between lateral static girder displacement and lateral static wind velocity; thus time series of lateral static girder displacement are decomposed into nonlinear correlation term and discreteness term, nonlinear correlation term of which is mathematically modeled bynth-order Fourier series with intervention of lateral static wind velocity and discreteness term of which is mathematically modeled by the combined models of ARMA(
p
,
q
) and EGARCH(
m
,
n
). Additionally, stable power spectrum density exists in time series of lateral dynamic girder displacement; thus time series of lateral dynamic girder displacement are mathematically modeled by harmonic superposition function. By comparison and verification between simulative and monitoring lateral displacements from September 1 to September 3, mathematical models are feasible and effective to simulate time series of lateral girder displacement from main girder of Sutong Cable-Stayed Bridge.


## 2. Bridge Monitoring and Sample Analysis

The bridge monitoring object for this research is the worldwide famous Sutong Cable-Stayed Bridge (in Jiangsu Province, China). Its whole structure form is single-spanned and double-hinged with the main span reaching 1088 m as shown in Figure1, and the main girder employs flat steel box type with 36.3 m wide and 4.0 m high as shown in Figure 2. 3D ultrasonic anemometers and GPS monitoring station are installed on two flanks of midspan cross-section from main girder (resp., shown in Figures 1 and 2) to continuously acquire wind data and displacement data with sample frequency of 1 Hz.Figure 1
Longitudinal layout of monitoring equipment on the Sutong Bridge (unit: m).Figure 2
Transverse layout of monitoring equipment on the flat steel box girder (unit: mm). (1) GPSMS: GPS monitoring station; (2) 3DUA: 3D ultrasonic anemometer.Specifically, the wind data from 3D ultrasonic anemometers embrace such three types as wind velocity, horizontal angle, and vertical angle in local coordinate system (Figure3), and the girder displacement data from GPS monitoring station contains absolute locations in WGS-84 coordinate system (Figure 3), which are supposed to deduct reference locations for analysis. Until now, the storage amount of monitoring data has increased to 93 million for each measurement point. Such considerable monitoring data cannot be totally applied to actual analysis; thus the monitoring data from upstream flank in the year 2012 are specially chosen.Figure 3
Description of the defined local and WGS-84 coordinate system. Notes: (1) local coordinate system:x
-
y
-
z; (2) WGS-84 coordinate system: X
-
Y
-
Z; (3) 90° of wind horizontal angle denotes the axis angle from x to y (within range of [0° 360°]); (4) 90° of wind vertical angle denotes axis angle from y to z (within range of [−90° 90°]); (5) 10.6° axis angle exists between axes Y and y of across the bridge.Taking the monitoring data in the whole August, for example, considering that lateral displacement effect of main girder is primarily reflected by wind load across the Sutong Bridge, hence time series of wind velocity and girder displacement are decomposed into they-axis in local coordinate system as shown in Figures 4(a) and 5(a), respectively, which can be obviously observed that either of two time series contain static variant trend in whole and dynamic stochastic fluctuation in part. Such two kinds of variation characteristics can be furthermore separated by 10-minute average process as shown in Figures 4(b), 4(c), 5(b), and 5(c), respectively.Time series of wind velocity alongy-axis in the whole August.
(a)
Time series of original wind velocity
(b)
Time series of static wind velocity
(c)
Time series of dynamic wind velocityTime series of girder displacements alongy-axis in the whole August.
(a)
Time series of girder displacement
(b)
Time series of static girder displacement
(c)
Time series of dynamic girder displacementBy comparison between Figures4(b) and 5(b), similar variation characteristics exist between static wind velocity and static girder displacement, which can be visually described by correlation scatter plots as shown in Figure 6, indicating apparent nonlinear correlation similar to quadratic parabolic curve. Therefore, static girder displacement can be mathematically expressed by static wind velocity, with consideration of definite discreteness affected by other random factors. Moreover, time series of dynamic girder displacement depict obvious steady stochastic fluctuation, as well as no variation of its power spectrum densities by time shown in Figure 7; thus dynamic girder displacement can be mathematically simulated by harmonic superposition method.Figure 6
Correlation scatter plots between static wind velocity and static girder displacement.Figure 7
Power spectrum densities of dynamic girder displacement in different periods.


## 3. Modeling Theory and Procedure

### 3.1. Modeling for Static Girder Displacement

#### 3.1.1. Modeling Theory and Method

Correlation scatter plots in Figure6 can be acquired by combination of nonlinear correlation term and discreteness term. Furthermore, nonlinear correlation term can be mathematically modeled by the n-order Fourier series; that is,
(1)
u
1
(
t
)
=
a
0
+
∑
k
=
1
n
(
a
k
cos
⁡
⁡
(
k
ω
v
(
t
)
)
+
b
k
sin
⁡
(
k
ω
v
(
t
)
)
)
,
where u
1
(
t
) denotes time series of fitting static displacement for nonlinear correlation term, v
(
t
) denotes time series of static wind velocity (Figure 4(b)), and a
k
,
b
k
(
k
=
0
,
1
,
…
,
n
), and ω are fitting parameters of nth-order Fourier series.With denotation of original static girder displacement asu
(
t
) (Figure 5(b)), u
(
t
) minus u
1
(
t
) acquires time series of discreteness u
2
(
t
), which contains three types of stochastic characteristics (autoregression, moving average, and heteroscedasticity) and can be mathematically described by the combined models of ARMA(
p
,
q
) and EGARCH(
m
,
n
). In detail, the ARMA(
p
,
q
) model defines the stochastic characteristics of autoregression and moving average as follows:(2)
u
2
(
t
)
=
c
+
∑
j
=
1
p
b
j
u
2
(
t
-
j
)
+
∑
j
=
0
q
c
j
ε
t
-
j
,
where c is the constant term, p and q, respectively, denote the orders of autoregression or moving average of u
2
(
t
), b
j and c
j, respectively, denote the coefficients of autoregression or moving average of u
2
(
t
) with c
0
=
1, and ε
t
-
j denotes the innovations process with time delay of j. Meanwhile, the other EGARCH(
m
,
n
) model defines the stochastic characteristic of heteroscedasticity as follows [9, 10]:
(3)
log
⁡
σ
t
2
=
κ
+
∑
k
=
1
m
g
k
log
⁡
σ
t
-
k
2
+
∑
k
=
1
n
v
k
[
|
ε
t
-
k
|
σ
t
-
k
-
E
{
|
z
t
-
k
|
}
]
+
∑
k
=
1
n
l
k
(
ε
t
-
k
σ
t
-
k
)
,
where σ
t denotes the conditional variance of the innovations process ε
t
,
κ is the constant term, m and n, respectively, denote the orders of the EGARCH(
m
,
n
) model, g
k
,
v
k, and l
k, respectively, denote the coefficients of the EGARCH(
m
,
n
) model, z
t is a standard, independent, and identically distributed random draw from some specified probability distribution such as Gaussian or Student’s t, and
(4)
E
{
|
z
t
-
k
|
}
=
E
(
|
ε
t
-
k
|
σ
t
-
k
)
=
{
2
π
Gaussian
ν
-
2
π
·
Γ
(
(
ν
-
1
)
/
2
)
Γ
(
ν
/
2
)
Student's
t
with degrees of freedom v
>
2. Considering that the EGARCH(
m
,
n
) model is treated as ARMA(
p
,
q
) models for log
⁡
σ
t
2, thus the stationarity constraint for the EGARCH(
m
,
n
) model is included by ensuring that the eigenvalues of the characteristic polynomial,
(5)
λ
m
-
G
1
λ
m
-
1
-
G
2
λ
m
-
2
-
⋯
-
G
m
,
are inside the unit circle, where λ and G
j are the variable and coefficients of the characteristic polynomial, respectively.During modeling process for time series of discretenessu
2
(
t
), the autocorrelation function ρ
(
s
) and partial correlation function β
(
s
) with their lag phase s are introduced for stationary test of u
2
(
t
). Specifically, the autocorrelation function ρ
(
s
) can be calculated as follows [11]:
(6)
ρ
(
s
)
=
(
1
/
N
)
∑
t
=
1
N
-
s
(
u
2
(
t
)
-
u
-
2
)
(
u
2
(
t
+
s
)
-
u
-
2
)
σ
ρ
2
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
s
=
0,1
,
2
,
…
,
N
-
1
,
where
(7)
u
-
2
=
1
N
∑
t
=
1
N
u
2
(
t
)
,
σ
ρ
2
=
1
N
∑
t
=
1
N
(
u
2
(
t
)
-
u
-
2
)
2
with N denoting the amount of u
2
(
t
). And the other partial correlation function β
(
s
) can be calculated through fitting successive autoregressive models of orders by ordinary least squares, retaining the last coefficient of each regression [12]. Besides, the AIC and BIC delimitation criteria are applied to determine model orders, with their statistical values a
i
c and b
i
c being, respectively, calculated as follows [13, 14]:
(8)
a
i
c
=
-
2
ln
⁡
(
L
)
+
2
K
1
,
b
i
c
=
-
2
ln
⁡
(
L
)
+
K
1
·
ln
⁡
(
K
2
)
,
where L denotes the optimized log-likelihood objective function (LLF) values associated with parameter estimates of the combined models, K
1 denotes the number of estimated parameters associated with each value in LLF, and K
2 denotes the sample size of the observed u
2
(
t
) associated with each LLF value.


#### 3.1.2. Detailed Procedure

Based on the theory and method above, the detailed procedure for mathematically modeling time series of static girder displacement is illustrated, taking the correlation scatter plots in the whole August in Figure6, for example, as follows.Step 1.
Fitting Fourier series for nonlinear correlation term. By means of the MATLAB fitting tools (utilizing the third-order Fourier series (1) to fit the correlation scatter plots) [10], the mathematical model of nonlinear correlation term is straightforward and acquired as shown in Figure 8(a), together with the estimated values of Fourier parameters presented in Table 1. By substitution of estimated values together with v
(
t
) into formula (1), time series of fitting static displacement u
1
(
t
) in the whole August are acquired as shown in Figure 8(b).Table 1
Estimated values of fitting parameters (-
12
<
v
(
t
)
<
21).
Fitting parameters
a
0
a
1
b
1
a
2
b
2
a
3
b
3
w
Estimated values
0.1102
−0.1026
0.0218
0.01214
0.006845
−0.00618
−0.003763
0.1374The fitting Fourier series and time series of fitting static displacement.
(a)
The fitting Fourier series of nonlinear correlation term
(b)
Time series of fitting static displacementStep 2.
Stationary test for time series of discreteness. Time series of discretenessu
2
(
t
) can be acquired by u
(
t
) minus u
1
(
t
) as shown in Figure 9(a). Autocorrelation function ρ
(
s
) and partial correlation function β
(
s
) of u
2
(
t
) are calculated with 50 lag phases shown in Figures 9(b) and 9(c), respectively, presenting that ρ
(
s
) is slowly converging into the 95% confidence intervals as the lag phase s increases, which indicates bad stationarity of u
2
(
t
) for mathematical modeling. Due to this, process of first-order difference for u
2
(
t
) is carried out as shown in Figure 10(a), with its ρ
(
s
) and β
(
s
) shown in Figures 10(b) and 10(c), both presenting rapid convergence into the 95% confidence intervals and verifying good stationarity for processed discreteness.Time series of original discreteness with itsρ
(
s
)
,
β
(
s
).
(a)
Time series of original discreteness
(b)
The autocorrelation functionρ
(
s
)
(c)
The partial correlation functionβ
(
s
)Time series of processed discreteness with itsρ
(
s
)
,
β
(
s
).
(a)
Time series of processed discreteness
(b)
The autocorrelation functionρ
(
s
)
(c)
The partial correlation functionβ
(
s
)Step 3.
Order determination of ARMA(
p
,
q
) and EGARCH(
m
,
n
). The orders of p and q are relative to the convergent forms of ρ
(
s
) and β
(
s
). That is, the ρ
(
s
) and β
(
s
) clearly present the trailing property in Figures 10(b) and 10(c), with 4th and 7th of lag phase s initially converging into the 95% confidence intervals, inferring that 7 and 4 are appropriately assigned to the orders of p and q, respectively [11, 12]. Furthermore, the orders of m and n are determined utilizing the AIC and BIC delimitation criterion. In detail, the statistical values of AIC and BIC from time series of processed discreteness (Figure 10(a)) are calculated, respectively, under integer assignment of m and n between 0 and 8 as shown in Figure 11, presenting that the most suitable orders of m and n are 2 and 1, respectively, corresponding to the minimum statistical values [13, 14].The statistical values of AIC and BIC
(a)
The statistical values of AIC
(b)
The statistical values of BICStep 4.
Parameter estimation and residual test of ordered models. Based on the specified orders above, model parameters (b
j, c
j, g
k, v
k, c, and κ) are further estimated to fit time series of processed discreteness, as shown in Tables 2 and 3, respectively. For testing the fitting effectiveness of estimated parameters, residuals between processed discreteness (Figure 10(a)) and its defined models with estimated parameters are analyzed using ρ
(
s
) and β
(
s
), as shown in Figure 12, which presents that both ρ
(
s
) and β
(
s
) are consistent in the 95% confidence intervals (except s is 0) and thus verifies good availability of estimated parameters for ordered models. Moreover, the standard errors between processed discreteness (Figure 10(a)) and its models with lower orders are shown in Figure 13.Table 2
Estimating model parameters of the ARMA(7, 4) model.
The combined models
The ARMA(7, 4) model
The fitting parameters
c
b
1
b
2
b
3
b
4
b
5
b
6
b
7
c
1
c
2
c
3
c
4
The values
0.00
−0.49
0.94
0.75
−0.32
−0.01
0.04
0.02
0.06
−1.26
−0.48
0.68Table 3
Estimating model parameters of the EGARCH(2, 1) model.
The combined models
The EGARCH(2, 1) model
The fitting parameters
κ
g
1
g
2
v
1
l
1
The values
−0.76
0.27
0.65
0.28
0.06The correlation functionsρ
(
s
) and β
(
s
) of residuals.
(a)
The autocorrelation functionρ
(
s
)
(b)
The partial correlation functionβ
(
s
)The standard errors between processed discreteness and its models with lower orders.
(a)
The standard errors whenm
=
1
,
n
=
1
(b)
The standard errors whenm
=
2
,
n
=
1Step 5.
Simulation for time series of discretenessu
2
(
t
) and static displacement u
(
t
). Based on the mathematical models of ARMA(
7,4
) and EGARCH(
2,1
) with estimated parameters, time series of simulative processed discreteness are shown in Figure 14(a), and through inverse calculation of first-order difference, time series of simulative original discreteness u
2
(
t
) are obtained as shown in Figure 14(b). Together with time series of fitting static displacement u
1
(
t
) in Figure 8(b), time series of simulative static displacement u
(
t
) are ultimately shown in Figure 14(c), definitely similar to time series of monitoring static displacement in Figure 5(b).Simulation for time series of discretenessu
2
(
t
) and static displacement u
(
t
).
(a)
Time series of simulative processed discreteness
(b)
Time series of simulative original discretenessu
2
(
t
)
(c)
Time series of simulative static displacementu
(
t
)


### 3.2. Modeling for Dynamic Girder Displacement

#### 3.2.1. Modeling Theory and Method

Time series of dynamic girder displacement show consistent power spectrum density (Figure7), which can be mathematically simulated by harmonic superposition method [15–17]. Primarily, function expression of power spectrum density should be confirmed to lay foundation for harmonic superposition. Considering that straightforward function fitting will ignore local feather of acute peak at 10−1 Hz around of frequency, thus the power spectrum density are decomposed into two parts: part one to fit the whole trend with ignorance of acute peak; part two to specifically fit the acute peak. Furthermore, each part can be expressed by the n-order Gaussian series in logarithmic form; that is,
(9)
lg
p
1
(
f
)
=
∑
k
=
1
n
a
k
e
-
(
(
lg
f
-
b
k
)
/
c
k
)
2
,
lg
p
2
(
f
)
=
∑
k
=
1
n
d
k
e
-
(
(
lg
f
-
r
k
)
/
h
k
)
2
,
where p
1
(
f
) denotes function expression of part one, p
2
(
f
) denotes function expression of part two, and a
k
,
b
k
,
c
k
,
d
k
,
h
k, and r
k are fitting parameters of p
1
(
f
) and p
2
(
f
)
,
k
=
0
,
1
,
…
,
n. Accordingly, the function expression of power spectrum density p
(
f
) can be formulized as follows:
(10)
p
(
f
)
=
p
1
(
f
)
+
p
2
(
f
)
.Based on power spectrum densityp
(
f
) above, variance σ
2 within frequency band of [
f
1
f
2
] can be calculated, utilizing the expanding characteristics of average power spectrum for stationary and stochastic process, as follows:
(11)
σ
2
=
∫
f
1
f
2
p
(
f
)
d
f
.
Then the frequency band of [f
1
f
2] is divided into n subintervals, and the whole values p
(
f
) in the ith subinterval can be represented by the corresponding value p
(
f
mid
,
i
) at the middle frequency f
mid
,
i; thus formula (11) is furthermore discretized as follows:
(12)
σ
2
≈
∑
i
=
1
n
p
(
f
mid
,
i
)
·
Δ
f
,
where
(13)
Δ
f
=
f
2
-
f
1
n
.
As for the ith subinterval, its sinusoidal function w
(
t
,
i
) containing middle frequency f
mid
,
i and standard deviation p
(
f
mid
,
i
)
·
Δ
f can be expressed as follows:
(14)
w
(
t
,
i
)
=
2
p
(
f
mid
,
i
)
·
Δ
f
sin
(
2
π
f
mid
,
i
t
+
θ
i
)
,
where θ
i is the random variable of uniform distribution within [
0
,
2
π
]. With superposition of w
(
t
,
i
) from each subinterval, the time series of dynamic girder displacement w
(
t
) can be mathematically modeled by harmonic superposition function as follows:
(15)
w
(
t
)
=
∑
i
=
1
n
2
p
(
f
mid
,
i
)
·
Δ
f
sin
(
2
π
f
mid
,
i
t
+
θ
i
)
.


#### 3.2.2. Detailed Procedure

Based on the theory and method above, the detailed procedure for mathematically modeled time series of dynamic girder displacement is illustrated, taking the correlation scatter plots from August 1 to August 8 in Figure15(a), for example, as follows.Original and processed power spectrum densities.
(a)
Original power spectrum density
(b)
Processed power spectrum densityStep 1.
Decreasing discreteness for original power spectrum density. Considering that certain discreteness existing in original power spectrum density can cover the acute peak at 10−1 Hz, adverse to fitting the n-order Gaussian series, therefore, average process with double frequency scales is carried out to decrease discreteness, specifically by 10−3 Hz average process within frequency band [10−1.5 Hz, 10−0.5 Hz] and 10−2 Hz average process within other frequency bands. Process result is shown in Figure 15(b), showing distinct acute peak for fitting the n-order Gaussian series.Step 2.
Fitting Gaussian series for two parts of processed power spectrum density. The processed power spectrum density can be divided into two parts: the whole trend with ignorance of acute peak as shown in Figure16(a); the specific acute peak as shown in Figure 16(b). By means of the MATLAB fitting tools (utilizing the fourth-order Gaussian series (9) to fit two parts) [10], the mathematical models p
1
(
f
) and p
2
(
f
) of two parts are, respectively, shown in Figures 15(a) and 15(b), together with estimated parameters presented in Table 4. The combination p
(
f
) of p
1
(
f
) and p
2
(
f
) is shown in Figure 16(c).Table 4
Estimated parameter values of the fourth-order Gaussian series.
Fitting curves
Estimated parameter values of the fourth-order Gaussian Series
a
1
a
2
a
3
a
4
b
1
b
2
b
3
b
4
c
1
c
2
c
3
c
4
p
1
(
f
)
−5.36
0.45
−4.86
−0.73
−4.98
−2.98
−0.75
−0.22
1.69
0.38
1.73
0.36
p
2
(
f
)
0.58
0.85
0.14
0.70
−0.96
−0.96
−1.20
−1.10
0.01
0.05
0.04
0.15Two parts of processed power spectrum densities and their fitting curves.
(a)
The first part and its fitting curvep
1
(
f
)
(b)
The second part and its fitting curvep
2
(
f
)
(c)
The combinationp
(
f
) of p
1
(
f
) and p
2
(
f
)Step 3.
Harmonic superposition for fitting power spectrum densityp
(
f
)
·
p
(
f
) is divided into 25000 subintervals within frequency bands [10−5 Hz, 10−0.5 Hz]. In each subinterval, one middle frequency f
mid
,
i and its corresponding value p
(
f
mid
,
i
) are existent and then substituted into (15) for summation; thus time series of dynamic girder displacement w
(
t
) can be acquired as shown in Figure 17(a) (simulated for 86400 s). By comparing its power spectrum density shown in Figure 17(b) with the monitoring one shown in Figure 15(a), good similarity of the whole trend and the acute peak verifies appropriateness of mathematical models for dynamic girder displacement.Time series of simulative dynamic displacementw
(
t
) and its power spectrum density.
(a)
Time series of simulative dynamic displacementw
(
t
)
(b)
Power spectrum density ofw
(
t
)


## 3.1. Modeling for Static Girder Displacement

### 3.1.1. Modeling Theory and Method

Correlation scatter plots in Figure6 can be acquired by combination of nonlinear correlation term and discreteness term. Furthermore, nonlinear correlation term can be mathematically modeled by the n-order Fourier series; that is,
(1)
u
1
(
t
)
=
a
0
+
∑
k
=
1
n
(
a
k
cos
⁡
⁡
(
k
ω
v
(
t
)
)
+
b
k
sin
⁡
(
k
ω
v
(
t
)
)
)
,
where u
1
(
t
) denotes time series of fitting static displacement for nonlinear correlation term, v
(
t
) denotes time series of static wind velocity (Figure 4(b)), and a
k
,
b
k
(
k
=
0
,
1
,
…
,
n
), and ω are fitting parameters of nth-order Fourier series.With denotation of original static girder displacement asu
(
t
) (Figure 5(b)), u
(
t
) minus u
1
(
t
) acquires time series of discreteness u
2
(
t
), which contains three types of stochastic characteristics (autoregression, moving average, and heteroscedasticity) and can be mathematically described by the combined models of ARMA(
p
,
q
) and EGARCH(
m
,
n
). In detail, the ARMA(
p
,
q
) model defines the stochastic characteristics of autoregression and moving average as follows:(2)
u
2
(
t
)
=
c
+
∑
j
=
1
p
b
j
u
2
(
t
-
j
)
+
∑
j
=
0
q
c
j
ε
t
-
j
,
where c is the constant term, p and q, respectively, denote the orders of autoregression or moving average of u
2
(
t
), b
j and c
j, respectively, denote the coefficients of autoregression or moving average of u
2
(
t
) with c
0
=
1, and ε
t
-
j denotes the innovations process with time delay of j. Meanwhile, the other EGARCH(
m
,
n
) model defines the stochastic characteristic of heteroscedasticity as follows [9, 10]:
(3)
log
⁡
σ
t
2
=
κ
+
∑
k
=
1
m
g
k
log
⁡
σ
t
-
k
2
+
∑
k
=
1
n
v
k
[
|
ε
t
-
k
|
σ
t
-
k
-
E
{
|
z
t
-
k
|
}
]
+
∑
k
=
1
n
l
k
(
ε
t
-
k
σ
t
-
k
)
,
where σ
t denotes the conditional variance of the innovations process ε
t
,
κ is the constant term, m and n, respectively, denote the orders of the EGARCH(
m
,
n
) model, g
k
,
v
k, and l
k, respectively, denote the coefficients of the EGARCH(
m
,
n
) model, z
t is a standard, independent, and identically distributed random draw from some specified probability distribution such as Gaussian or Student’s t, and
(4)
E
{
|
z
t
-
k
|
}
=
E
(
|
ε
t
-
k
|
σ
t
-
k
)
=
{
2
π
Gaussian
ν
-
2
π
·
Γ
(
(
ν
-
1
)
/
2
)
Γ
(
ν
/
2
)
Student's
t
with degrees of freedom v
>
2. Considering that the EGARCH(
m
,
n
) model is treated as ARMA(
p
,
q
) models for log
⁡
σ
t
2, thus the stationarity constraint for the EGARCH(
m
,
n
) model is included by ensuring that the eigenvalues of the characteristic polynomial,
(5)
λ
m
-
G
1
λ
m
-
1
-
G
2
λ
m
-
2
-
⋯
-
G
m
,
are inside the unit circle, where λ and G
j are the variable and coefficients of the characteristic polynomial, respectively.During modeling process for time series of discretenessu
2
(
t
), the autocorrelation function ρ
(
s
) and partial correlation function β
(
s
) with their lag phase s are introduced for stationary test of u
2
(
t
). Specifically, the autocorrelation function ρ
(
s
) can be calculated as follows [11]:
(6)
ρ
(
s
)
=
(
1
/
N
)
∑
t
=
1
N
-
s
(
u
2
(
t
)
-
u
-
2
)
(
u
2
(
t
+
s
)
-
u
-
2
)
σ
ρ
2
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
s
=
0,1
,
2
,
…
,
N
-
1
,
where
(7)
u
-
2
=
1
N
∑
t
=
1
N
u
2
(
t
)
,
σ
ρ
2
=
1
N
∑
t
=
1
N
(
u
2
(
t
)
-
u
-
2
)
2
with N denoting the amount of u
2
(
t
). And the other partial correlation function β
(
s
) can be calculated through fitting successive autoregressive models of orders by ordinary least squares, retaining the last coefficient of each regression [12]. Besides, the AIC and BIC delimitation criteria are applied to determine model orders, with their statistical values a
i
c and b
i
c being, respectively, calculated as follows [13, 14]:
(8)
a
i
c
=
-
2
ln
⁡
(
L
)
+
2
K
1
,
b
i
c
=
-
2
ln
⁡
(
L
)
+
K
1
·
ln
⁡
(
K
2
)
,
where L denotes the optimized log-likelihood objective function (LLF) values associated with parameter estimates of the combined models, K
1 denotes the number of estimated parameters associated with each value in LLF, and K
2 denotes the sample size of the observed u
2
(
t
) associated with each LLF value.


### 3.1.2. Detailed Procedure

Based on the theory and method above, the detailed procedure for mathematically modeling time series of static girder displacement is illustrated, taking the correlation scatter plots in the whole August in Figure6, for example, as follows.Step 1.
Fitting Fourier series for nonlinear correlation term. By means of the MATLAB fitting tools (utilizing the third-order Fourier series (1) to fit the correlation scatter plots) [10], the mathematical model of nonlinear correlation term is straightforward and acquired as shown in Figure 8(a), together with the estimated values of Fourier parameters presented in Table 1. By substitution of estimated values together with v
(
t
) into formula (1), time series of fitting static displacement u
1
(
t
) in the whole August are acquired as shown in Figure 8(b).Table 1
Estimated values of fitting parameters (-
12
<
v
(
t
)
<
21).
Fitting parameters
a
0
a
1
b
1
a
2
b
2
a
3
b
3
w
Estimated values
0.1102
−0.1026
0.0218
0.01214
0.006845
−0.00618
−0.003763
0.1374The fitting Fourier series and time series of fitting static displacement.
(a)
The fitting Fourier series of nonlinear correlation term
(b)
Time series of fitting static displacementStep 2.
Stationary test for time series of discreteness. Time series of discretenessu
2
(
t
) can be acquired by u
(
t
) minus u
1
(
t
) as shown in Figure 9(a). Autocorrelation function ρ
(
s
) and partial correlation function β
(
s
) of u
2
(
t
) are calculated with 50 lag phases shown in Figures 9(b) and 9(c), respectively, presenting that ρ
(
s
) is slowly converging into the 95% confidence intervals as the lag phase s increases, which indicates bad stationarity of u
2
(
t
) for mathematical modeling. Due to this, process of first-order difference for u
2
(
t
) is carried out as shown in Figure 10(a), with its ρ
(
s
) and β
(
s
) shown in Figures 10(b) and 10(c), both presenting rapid convergence into the 95% confidence intervals and verifying good stationarity for processed discreteness.Time series of original discreteness with itsρ
(
s
)
,
β
(
s
).
(a)
Time series of original discreteness
(b)
The autocorrelation functionρ
(
s
)
(c)
The partial correlation functionβ
(
s
)Time series of processed discreteness with itsρ
(
s
)
,
β
(
s
).
(a)
Time series of processed discreteness
(b)
The autocorrelation functionρ
(
s
)
(c)
The partial correlation functionβ
(
s
)Step 3.
Order determination of ARMA(
p
,
q
) and EGARCH(
m
,
n
). The orders of p and q are relative to the convergent forms of ρ
(
s
) and β
(
s
). That is, the ρ
(
s
) and β
(
s
) clearly present the trailing property in Figures 10(b) and 10(c), with 4th and 7th of lag phase s initially converging into the 95% confidence intervals, inferring that 7 and 4 are appropriately assigned to the orders of p and q, respectively [11, 12]. Furthermore, the orders of m and n are determined utilizing the AIC and BIC delimitation criterion. In detail, the statistical values of AIC and BIC from time series of processed discreteness (Figure 10(a)) are calculated, respectively, under integer assignment of m and n between 0 and 8 as shown in Figure 11, presenting that the most suitable orders of m and n are 2 and 1, respectively, corresponding to the minimum statistical values [13, 14].The statistical values of AIC and BIC
(a)
The statistical values of AIC
(b)
The statistical values of BICStep 4.
Parameter estimation and residual test of ordered models. Based on the specified orders above, model parameters (b
j, c
j, g
k, v
k, c, and κ) are further estimated to fit time series of processed discreteness, as shown in Tables 2 and 3, respectively. For testing the fitting effectiveness of estimated parameters, residuals between processed discreteness (Figure 10(a)) and its defined models with estimated parameters are analyzed using ρ
(
s
) and β
(
s
), as shown in Figure 12, which presents that both ρ
(
s
) and β
(
s
) are consistent in the 95% confidence intervals (except s is 0) and thus verifies good availability of estimated parameters for ordered models. Moreover, the standard errors between processed discreteness (Figure 10(a)) and its models with lower orders are shown in Figure 13.Table 2
Estimating model parameters of the ARMA(7, 4) model.
The combined models
The ARMA(7, 4) model
The fitting parameters
c
b
1
b
2
b
3
b
4
b
5
b
6
b
7
c
1
c
2
c
3
c
4
The values
0.00
−0.49
0.94
0.75
−0.32
−0.01
0.04
0.02
0.06
−1.26
−0.48
0.68Table 3
Estimating model parameters of the EGARCH(2, 1) model.
The combined models
The EGARCH(2, 1) model
The fitting parameters
κ
g
1
g
2
v
1
l
1
The values
−0.76
0.27
0.65
0.28
0.06The correlation functionsρ
(
s
) and β
(
s
) of residuals.
(a)
The autocorrelation functionρ
(
s
)
(b)
The partial correlation functionβ
(
s
)The standard errors between processed discreteness and its models with lower orders.
(a)
The standard errors whenm
=
1
,
n
=
1
(b)
The standard errors whenm
=
2
,
n
=
1Step 5.
Simulation for time series of discretenessu
2
(
t
) and static displacement u
(
t
). Based on the mathematical models of ARMA(
7,4
) and EGARCH(
2,1
) with estimated parameters, time series of simulative processed discreteness are shown in Figure 14(a), and through inverse calculation of first-order difference, time series of simulative original discreteness u
2
(
t
) are obtained as shown in Figure 14(b). Together with time series of fitting static displacement u
1
(
t
) in Figure 8(b), time series of simulative static displacement u
(
t
) are ultimately shown in Figure 14(c), definitely similar to time series of monitoring static displacement in Figure 5(b).Simulation for time series of discretenessu
2
(
t
) and static displacement u
(
t
).
(a)
Time series of simulative processed discreteness
(b)
Time series of simulative original discretenessu
2
(
t
)
(c)
Time series of simulative static displacementu
(
t
)


## 3.1.1. Modeling Theory and Method

Correlation scatter plots in Figure6 can be acquired by combination of nonlinear correlation term and discreteness term. Furthermore, nonlinear correlation term can be mathematically modeled by the n-order Fourier series; that is,
(1)
u
1
(
t
)
=
a
0
+
∑
k
=
1
n
(
a
k
cos
⁡
⁡
(
k
ω
v
(
t
)
)
+
b
k
sin
⁡
(
k
ω
v
(
t
)
)
)
,
where u
1
(
t
) denotes time series of fitting static displacement for nonlinear correlation term, v
(
t
) denotes time series of static wind velocity (Figure 4(b)), and a
k
,
b
k
(
k
=
0
,
1
,
…
,
n
), and ω are fitting parameters of nth-order Fourier series.With denotation of original static girder displacement asu
(
t
) (Figure 5(b)), u
(
t
) minus u
1
(
t
) acquires time series of discreteness u
2
(
t
), which contains three types of stochastic characteristics (autoregression, moving average, and heteroscedasticity) and can be mathematically described by the combined models of ARMA(
p
,
q
) and EGARCH(
m
,
n
). In detail, the ARMA(
p
,
q
) model defines the stochastic characteristics of autoregression and moving average as follows:(2)
u
2
(
t
)
=
c
+
∑
j
=
1
p
b
j
u
2
(
t
-
j
)
+
∑
j
=
0
q
c
j
ε
t
-
j
,
where c is the constant term, p and q, respectively, denote the orders of autoregression or moving average of u
2
(
t
), b
j and c
j, respectively, denote the coefficients of autoregression or moving average of u
2
(
t
) with c
0
=
1, and ε
t
-
j denotes the innovations process with time delay of j. Meanwhile, the other EGARCH(
m
,
n
) model defines the stochastic characteristic of heteroscedasticity as follows [9, 10]:
(3)
log
⁡
σ
t
2
=
κ
+
∑
k
=
1
m
g
k
log
⁡
σ
t
-
k
2
+
∑
k
=
1
n
v
k
[
|
ε
t
-
k
|
σ
t
-
k
-
E
{
|
z
t
-
k
|
}
]
+
∑
k
=
1
n
l
k
(
ε
t
-
k
σ
t
-
k
)
,
where σ
t denotes the conditional variance of the innovations process ε
t
,
κ is the constant term, m and n, respectively, denote the orders of the EGARCH(
m
,
n
) model, g
k
,
v
k, and l
k, respectively, denote the coefficients of the EGARCH(
m
,
n
) model, z
t is a standard, independent, and identically distributed random draw from some specified probability distribution such as Gaussian or Student’s t, and
(4)
E
{
|
z
t
-
k
|
}
=
E
(
|
ε
t
-
k
|
σ
t
-
k
)
=
{
2
π
Gaussian
ν
-
2
π
·
Γ
(
(
ν
-
1
)
/
2
)
Γ
(
ν
/
2
)
Student's
t
with degrees of freedom v
>
2. Considering that the EGARCH(
m
,
n
) model is treated as ARMA(
p
,
q
) models for log
⁡
σ
t
2, thus the stationarity constraint for the EGARCH(
m
,
n
) model is included by ensuring that the eigenvalues of the characteristic polynomial,
(5)
λ
m
-
G
1
λ
m
-
1
-
G
2
λ
m
-
2
-
⋯
-
G
m
,
are inside the unit circle, where λ and G
j are the variable and coefficients of the characteristic polynomial, respectively.During modeling process for time series of discretenessu
2
(
t
), the autocorrelation function ρ
(
s
) and partial correlation function β
(
s
) with their lag phase s are introduced for stationary test of u
2
(
t
). Specifically, the autocorrelation function ρ
(
s
) can be calculated as follows [11]:
(6)
ρ
(
s
)
=
(
1
/
N
)
∑
t
=
1
N
-
s
(
u
2
(
t
)
-
u
-
2
)
(
u
2
(
t
+
s
)
-
u
-
2
)
σ
ρ
2
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
A
s
=
0,1
,
2
,
…
,
N
-
1
,
where
(7)
u
-
2
=
1
N
∑
t
=
1
N
u
2
(
t
)
,
σ
ρ
2
=
1
N
∑
t
=
1
N
(
u
2
(
t
)
-
u
-
2
)
2
with N denoting the amount of u
2
(
t
). And the other partial correlation function β
(
s
) can be calculated through fitting successive autoregressive models of orders by ordinary least squares, retaining the last coefficient of each regression [12]. Besides, the AIC and BIC delimitation criteria are applied to determine model orders, with their statistical values a
i
c and b
i
c being, respectively, calculated as follows [13, 14]:
(8)
a
i
c
=
-
2
ln
⁡
(
L
)
+
2
K
1
,
b
i
c
=
-
2
ln
⁡
(
L
)
+
K
1
·
ln
⁡
(
K
2
)
,
where L denotes the optimized log-likelihood objective function (LLF) values associated with parameter estimates of the combined models, K
1 denotes the number of estimated parameters associated with each value in LLF, and K
2 denotes the sample size of the observed u
2
(
t
) associated with each LLF value.


## 3.1.2. Detailed Procedure

Based on the theory and method above, the detailed procedure for mathematically modeling time series of static girder displacement is illustrated, taking the correlation scatter plots in the whole August in Figure6, for example, as follows.Step 1.
Fitting Fourier series for nonlinear correlation term. By means of the MATLAB fitting tools (utilizing the third-order Fourier series (1) to fit the correlation scatter plots) [10], the mathematical model of nonlinear correlation term is straightforward and acquired as shown in Figure 8(a), together with the estimated values of Fourier parameters presented in Table 1. By substitution of estimated values together with v
(
t
) into formula (1), time series of fitting static displacement u
1
(
t
) in the whole August are acquired as shown in Figure 8(b).Table 1
Estimated values of fitting parameters (-
12
<
v
(
t
)
<
21).
Fitting parameters
a
0
a
1
b
1
a
2
b
2
a
3
b
3
w
Estimated values
0.1102
−0.1026
0.0218
0.01214
0.006845
−0.00618
−0.003763
0.1374The fitting Fourier series and time series of fitting static displacement.
(a)
The fitting Fourier series of nonlinear correlation term
(b)
Time series of fitting static displacementStep 2.
Stationary test for time series of discreteness. Time series of discretenessu
2
(
t
) can be acquired by u
(
t
) minus u
1
(
t
) as shown in Figure 9(a). Autocorrelation function ρ
(
s
) and partial correlation function β
(
s
) of u
2
(
t
) are calculated with 50 lag phases shown in Figures 9(b) and 9(c), respectively, presenting that ρ
(
s
) is slowly converging into the 95% confidence intervals as the lag phase s increases, which indicates bad stationarity of u
2
(
t
) for mathematical modeling. Due to this, process of first-order difference for u
2
(
t
) is carried out as shown in Figure 10(a), with its ρ
(
s
) and β
(
s
) shown in Figures 10(b) and 10(c), both presenting rapid convergence into the 95% confidence intervals and verifying good stationarity for processed discreteness.Time series of original discreteness with itsρ
(
s
)
,
β
(
s
).
(a)
Time series of original discreteness
(b)
The autocorrelation functionρ
(
s
)
(c)
The partial correlation functionβ
(
s
)Time series of processed discreteness with itsρ
(
s
)
,
β
(
s
).
(a)
Time series of processed discreteness
(b)
The autocorrelation functionρ
(
s
)
(c)
The partial correlation functionβ
(
s
)Step 3.
Order determination of ARMA(
p
,
q
) and EGARCH(
m
,
n
). The orders of p and q are relative to the convergent forms of ρ
(
s
) and β
(
s
). That is, the ρ
(
s
) and β
(
s
) clearly present the trailing property in Figures 10(b) and 10(c), with 4th and 7th of lag phase s initially converging into the 95% confidence intervals, inferring that 7 and 4 are appropriately assigned to the orders of p and q, respectively [11, 12]. Furthermore, the orders of m and n are determined utilizing the AIC and BIC delimitation criterion. In detail, the statistical values of AIC and BIC from time series of processed discreteness (Figure 10(a)) are calculated, respectively, under integer assignment of m and n between 0 and 8 as shown in Figure 11, presenting that the most suitable orders of m and n are 2 and 1, respectively, corresponding to the minimum statistical values [13, 14].The statistical values of AIC and BIC
(a)
The statistical values of AIC
(b)
The statistical values of BICStep 4.
Parameter estimation and residual test of ordered models. Based on the specified orders above, model parameters (b
j, c
j, g
k, v
k, c, and κ) are further estimated to fit time series of processed discreteness, as shown in Tables 2 and 3, respectively. For testing the fitting effectiveness of estimated parameters, residuals between processed discreteness (Figure 10(a)) and its defined models with estimated parameters are analyzed using ρ
(
s
) and β
(
s
), as shown in Figure 12, which presents that both ρ
(
s
) and β
(
s
) are consistent in the 95% confidence intervals (except s is 0) and thus verifies good availability of estimated parameters for ordered models. Moreover, the standard errors between processed discreteness (Figure 10(a)) and its models with lower orders are shown in Figure 13.Table 2
Estimating model parameters of the ARMA(7, 4) model.
The combined models
The ARMA(7, 4) model
The fitting parameters
c
b
1
b
2
b
3
b
4
b
5
b
6
b
7
c
1
c
2
c
3
c
4
The values
0.00
−0.49
0.94
0.75
−0.32
−0.01
0.04
0.02
0.06
−1.26
−0.48
0.68Table 3
Estimating model parameters of the EGARCH(2, 1) model.
The combined models
The EGARCH(2, 1) model
The fitting parameters
κ
g
1
g
2
v
1
l
1
The values
−0.76
0.27
0.65
0.28
0.06The correlation functionsρ
(
s
) and β
(
s
) of residuals.
(a)
The autocorrelation functionρ
(
s
)
(b)
The partial correlation functionβ
(
s
)The standard errors between processed discreteness and its models with lower orders.
(a)
The standard errors whenm
=
1
,
n
=
1
(b)
The standard errors whenm
=
2
,
n
=
1Step 5.
Simulation for time series of discretenessu
2
(
t
) and static displacement u
(
t
). Based on the mathematical models of ARMA(
7,4
) and EGARCH(
2,1
) with estimated parameters, time series of simulative processed discreteness are shown in Figure 14(a), and through inverse calculation of first-order difference, time series of simulative original discreteness u
2
(
t
) are obtained as shown in Figure 14(b). Together with time series of fitting static displacement u
1
(
t
) in Figure 8(b), time series of simulative static displacement u
(
t
) are ultimately shown in Figure 14(c), definitely similar to time series of monitoring static displacement in Figure 5(b).Simulation for time series of discretenessu
2
(
t
) and static displacement u
(
t
).
(a)
Time series of simulative processed discreteness
(b)
Time series of simulative original discretenessu
2
(
t
)
(c)
Time series of simulative static displacementu
(
t
)


## 3.2. Modeling for Dynamic Girder Displacement

### 3.2.1. Modeling Theory and Method

Time series of dynamic girder displacement show consistent power spectrum density (Figure7), which can be mathematically simulated by harmonic superposition method [15–17]. Primarily, function expression of power spectrum density should be confirmed to lay foundation for harmonic superposition. Considering that straightforward function fitting will ignore local feather of acute peak at 10−1 Hz around of frequency, thus the power spectrum density are decomposed into two parts: part one to fit the whole trend with ignorance of acute peak; part two to specifically fit the acute peak. Furthermore, each part can be expressed by the n-order Gaussian series in logarithmic form; that is,
(9)
lg
p
1
(
f
)
=
∑
k
=
1
n
a
k
e
-
(
(
lg
f
-
b
k
)
/
c
k
)
2
,
lg
p
2
(
f
)
=
∑
k
=
1
n
d
k
e
-
(
(
lg
f
-
r
k
)
/
h
k
)
2
,
where p
1
(
f
) denotes function expression of part one, p
2
(
f
) denotes function expression of part two, and a
k
,
b
k
,
c
k
,
d
k
,
h
k, and r
k are fitting parameters of p
1
(
f
) and p
2
(
f
)
,
k
=
0
,
1
,
…
,
n. Accordingly, the function expression of power spectrum density p
(
f
) can be formulized as follows:
(10)
p
(
f
)
=
p
1
(
f
)
+
p
2
(
f
)
.Based on power spectrum densityp
(
f
) above, variance σ
2 within frequency band of [
f
1
f
2
] can be calculated, utilizing the expanding characteristics of average power spectrum for stationary and stochastic process, as follows:
(11)
σ
2
=
∫
f
1
f
2
p
(
f
)
d
f
.
Then the frequency band of [f
1
f
2] is divided into n subintervals, and the whole values p
(
f
) in the ith subinterval can be represented by the corresponding value p
(
f
mid
,
i
) at the middle frequency f
mid
,
i; thus formula (11) is furthermore discretized as follows:
(12)
σ
2
≈
∑
i
=
1
n
p
(
f
mid
,
i
)
·
Δ
f
,
where
(13)
Δ
f
=
f
2
-
f
1
n
.
As for the ith subinterval, its sinusoidal function w
(
t
,
i
) containing middle frequency f
mid
,
i and standard deviation p
(
f
mid
,
i
)
·
Δ
f can be expressed as follows:
(14)
w
(
t
,
i
)
=
2
p
(
f
mid
,
i
)
·
Δ
f
sin
(
2
π
f
mid
,
i
t
+
θ
i
)
,
where θ
i is the random variable of uniform distribution within [
0
,
2
π
]. With superposition of w
(
t
,
i
) from each subinterval, the time series of dynamic girder displacement w
(
t
) can be mathematically modeled by harmonic superposition function as follows:
(15)
w
(
t
)
=
∑
i
=
1
n
2
p
(
f
mid
,
i
)
·
Δ
f
sin
(
2
π
f
mid
,
i
t
+
θ
i
)
.


### 3.2.2. Detailed Procedure

Based on the theory and method above, the detailed procedure for mathematically modeled time series of dynamic girder displacement is illustrated, taking the correlation scatter plots from August 1 to August 8 in Figure15(a), for example, as follows.Original and processed power spectrum densities.
(a)
Original power spectrum density
(b)
Processed power spectrum densityStep 1.
Decreasing discreteness for original power spectrum density. Considering that certain discreteness existing in original power spectrum density can cover the acute peak at 10−1 Hz, adverse to fitting the n-order Gaussian series, therefore, average process with double frequency scales is carried out to decrease discreteness, specifically by 10−3 Hz average process within frequency band [10−1.5 Hz, 10−0.5 Hz] and 10−2 Hz average process within other frequency bands. Process result is shown in Figure 15(b), showing distinct acute peak for fitting the n-order Gaussian series.Step 2.
Fitting Gaussian series for two parts of processed power spectrum density. The processed power spectrum density can be divided into two parts: the whole trend with ignorance of acute peak as shown in Figure16(a); the specific acute peak as shown in Figure 16(b). By means of the MATLAB fitting tools (utilizing the fourth-order Gaussian series (9) to fit two parts) [10], the mathematical models p
1
(
f
) and p
2
(
f
) of two parts are, respectively, shown in Figures 15(a) and 15(b), together with estimated parameters presented in Table 4. The combination p
(
f
) of p
1
(
f
) and p
2
(
f
) is shown in Figure 16(c).Table 4
Estimated parameter values of the fourth-order Gaussian series.
Fitting curves
Estimated parameter values of the fourth-order Gaussian Series
a
1
a
2
a
3
a
4
b
1
b
2
b
3
b
4
c
1
c
2
c
3
c
4
p
1
(
f
)
−5.36
0.45
−4.86
−0.73
−4.98
−2.98
−0.75
−0.22
1.69
0.38
1.73
0.36
p
2
(
f
)
0.58
0.85
0.14
0.70
−0.96
−0.96
−1.20
−1.10
0.01
0.05
0.04
0.15Two parts of processed power spectrum densities and their fitting curves.
(a)
The first part and its fitting curvep
1
(
f
)
(b)
The second part and its fitting curvep
2
(
f
)
(c)
The combinationp
(
f
) of p
1
(
f
) and p
2
(
f
)Step 3.
Harmonic superposition for fitting power spectrum densityp
(
f
)
·
p
(
f
) is divided into 25000 subintervals within frequency bands [10−5 Hz, 10−0.5 Hz]. In each subinterval, one middle frequency f
mid
,
i and its corresponding value p
(
f
mid
,
i
) are existent and then substituted into (15) for summation; thus time series of dynamic girder displacement w
(
t
) can be acquired as shown in Figure 17(a) (simulated for 86400 s). By comparing its power spectrum density shown in Figure 17(b) with the monitoring one shown in Figure 15(a), good similarity of the whole trend and the acute peak verifies appropriateness of mathematical models for dynamic girder displacement.Time series of simulative dynamic displacementw
(
t
) and its power spectrum density.
(a)
Time series of simulative dynamic displacementw
(
t
)
(b)
Power spectrum density ofw
(
t
)


## 3.2.1. Modeling Theory and Method

Time series of dynamic girder displacement show consistent power spectrum density (Figure7), which can be mathematically simulated by harmonic superposition method [15–17]. Primarily, function expression of power spectrum density should be confirmed to lay foundation for harmonic superposition. Considering that straightforward function fitting will ignore local feather of acute peak at 10−1 Hz around of frequency, thus the power spectrum density are decomposed into two parts: part one to fit the whole trend with ignorance of acute peak; part two to specifically fit the acute peak. Furthermore, each part can be expressed by the n-order Gaussian series in logarithmic form; that is,
(9)
lg
p
1
(
f
)
=
∑
k
=
1
n
a
k
e
-
(
(
lg
f
-
b
k
)
/
c
k
)
2
,
lg
p
2
(
f
)
=
∑
k
=
1
n
d
k
e
-
(
(
lg
f
-
r
k
)
/
h
k
)
2
,
where p
1
(
f
) denotes function expression of part one, p
2
(
f
) denotes function expression of part two, and a
k
,
b
k
,
c
k
,
d
k
,
h
k, and r
k are fitting parameters of p
1
(
f
) and p
2
(
f
)
,
k
=
0
,
1
,
…
,
n. Accordingly, the function expression of power spectrum density p
(
f
) can be formulized as follows:
(10)
p
(
f
)
=
p
1
(
f
)
+
p
2
(
f
)
.Based on power spectrum densityp
(
f
) above, variance σ
2 within frequency band of [
f
1
f
2
] can be calculated, utilizing the expanding characteristics of average power spectrum for stationary and stochastic process, as follows:
(11)
σ
2
=
∫
f
1
f
2
p
(
f
)
d
f
.
Then the frequency band of [f
1
f
2] is divided into n subintervals, and the whole values p
(
f
) in the ith subinterval can be represented by the corresponding value p
(
f
mid
,
i
) at the middle frequency f
mid
,
i; thus formula (11) is furthermore discretized as follows:
(12)
σ
2
≈
∑
i
=
1
n
p
(
f
mid
,
i
)
·
Δ
f
,
where
(13)
Δ
f
=
f
2
-
f
1
n
.
As for the ith subinterval, its sinusoidal function w
(
t
,
i
) containing middle frequency f
mid
,
i and standard deviation p
(
f
mid
,
i
)
·
Δ
f can be expressed as follows:
(14)
w
(
t
,
i
)
=
2
p
(
f
mid
,
i
)
·
Δ
f
sin
(
2
π
f
mid
,
i
t
+
θ
i
)
,
where θ
i is the random variable of uniform distribution within [
0
,
2
π
]. With superposition of w
(
t
,
i
) from each subinterval, the time series of dynamic girder displacement w
(
t
) can be mathematically modeled by harmonic superposition function as follows:
(15)
w
(
t
)
=
∑
i
=
1
n
2
p
(
f
mid
,
i
)
·
Δ
f
sin
(
2
π
f
mid
,
i
t
+
θ
i
)
.


## 3.2.2. Detailed Procedure

Based on the theory and method above, the detailed procedure for mathematically modeled time series of dynamic girder displacement is illustrated, taking the correlation scatter plots from August 1 to August 8 in Figure15(a), for example, as follows.Original and processed power spectrum densities.
(a)
Original power spectrum density
(b)
Processed power spectrum densityStep 1.
Decreasing discreteness for original power spectrum density. Considering that certain discreteness existing in original power spectrum density can cover the acute peak at 10−1 Hz, adverse to fitting the n-order Gaussian series, therefore, average process with double frequency scales is carried out to decrease discreteness, specifically by 10−3 Hz average process within frequency band [10−1.5 Hz, 10−0.5 Hz] and 10−2 Hz average process within other frequency bands. Process result is shown in Figure 15(b), showing distinct acute peak for fitting the n-order Gaussian series.Step 2.
Fitting Gaussian series for two parts of processed power spectrum density. The processed power spectrum density can be divided into two parts: the whole trend with ignorance of acute peak as shown in Figure16(a); the specific acute peak as shown in Figure 16(b). By means of the MATLAB fitting tools (utilizing the fourth-order Gaussian series (9) to fit two parts) [10], the mathematical models p
1
(
f
) and p
2
(
f
) of two parts are, respectively, shown in Figures 15(a) and 15(b), together with estimated parameters presented in Table 4. The combination p
(
f
) of p
1
(
f
) and p
2
(
f
) is shown in Figure 16(c).Table 4
Estimated parameter values of the fourth-order Gaussian series.
Fitting curves
Estimated parameter values of the fourth-order Gaussian Series
a
1
a
2
a
3
a
4
b
1
b
2
b
3
b
4
c
1
c
2
c
3
c
4
p
1
(
f
)
−5.36
0.45
−4.86
−0.73
−4.98
−2.98
−0.75
−0.22
1.69
0.38
1.73
0.36
p
2
(
f
)
0.58
0.85
0.14
0.70
−0.96
−0.96
−1.20
−1.10
0.01
0.05
0.04
0.15Two parts of processed power spectrum densities and their fitting curves.
(a)
The first part and its fitting curvep
1
(
f
)
(b)
The second part and its fitting curvep
2
(
f
)
(c)
The combinationp
(
f
) of p
1
(
f
) and p
2
(
f
)Step 3.
Harmonic superposition for fitting power spectrum densityp
(
f
)
·
p
(
f
) is divided into 25000 subintervals within frequency bands [10−5 Hz, 10−0.5 Hz]. In each subinterval, one middle frequency f
mid
,
i and its corresponding value p
(
f
mid
,
i
) are existent and then substituted into (15) for summation; thus time series of dynamic girder displacement w
(
t
) can be acquired as shown in Figure 17(a) (simulated for 86400 s). By comparing its power spectrum density shown in Figure 17(b) with the monitoring one shown in Figure 15(a), good similarity of the whole trend and the acute peak verifies appropriateness of mathematical models for dynamic girder displacement.Time series of simulative dynamic displacementw
(
t
) and its power spectrum density.
(a)
Time series of simulative dynamic displacementw
(
t
)
(b)
Power spectrum density ofw
(
t
)


## 4. Model Test and Evaluation

According to mathematical modeling process above, time series of lateral girder displacement are expressed by combination of third-order Fourier seriesu
1
(
t
), ARMA(7,4), EGARCH(
2,1
), and harmonic superposition function w
(
t
). For verifying feasibility and effectiveness of the whole mathematical models, time series from September 1 to September 3 are simulated with intervention of monitoring static wind (Figure 18(a)) and then compared with monitoring ones during same period (Figure 18(b)).Time series of simulative and monitoring lateral girder displacement.
(a)
Time series of simulative lateral displacement
(b)
Time series of monitoring lateral displacementAccording to mathematical modeling theory and procedure above, comparison is divided into two parts: (1) time series of static girder displacement; (2) time series of dynamic girder displacement. As for the first part, its simulative and monitoring results are shown in Figure19(a) and linear fitting curve of correlation scatter plots is shown in Figure 19(b), presenting consistent variation tendency in Figure 19(a) and d
s
(
t
) approximating to d
m
(
t
) in Figure 19(b) (where d
s
(
t
) and d
m
(
t
), respectively, denote simulative and monitoring results), which verifies good feasibility and effectiveness of mathematical models for static girder displacement. As for the second part, power spectrum densities of simulative and monitoring results are shown in Figure 20, presenting uniform variation tendency of both whole trends and acute peaks, which verifies good feasibility and effectiveness of mathematical models for dynamic girder displacement. Therefore, mathematical models above can be reasonably utilized to simulate time series of lateral girder displacement from main girder of Sutong Cable-Stayed Bridge.Simulative and monitoring results and their correlation scatter plots with fitting curve.
(a)
Simulative and monitoring results
(b)
Correlation scatter plots with fitting curveFigure 20
Power spectrum densities of simulative and monitoring results.


## 5. Conclusions

Based on monitoring data from main girder of Sutong Cable-Stayed Bridge, time series of lateral girder displacement effect are mathematically modeled by methods of fitting Fourier series and Gaussian series, combined models of ARMA(7,4) and EGARCH(2,1), and harmonic superposition function. And conclusions can be drawn as follows.(1)
Scatter plots between lateral static wind velocity and lateral static displacement present apparent nonlinear correlation, which is similar to quadratic parabolic curve, and time series of lateral dynamic displacement contain obvious stable power spectrum density with no variation by time.
(2)
Time series of lateral static displacement can be decomposed into nonlinear correlation term and discreteness term. Moreover, nonlinear correlation term can be mathematically modeled by third-order Fourier series with intervention of lateral static wind velocity, and discreteness term can be mathematically modeled by the combined models of ARMA(
7,4
) and EGARCH(
2,1
).
(3)
Through decreasing discreteness in double frequency scales and division in double frequency bands, power spectrum density of lateral dynamic displacement can be mathematically modeled by the fourth-order Gaussian series, and time series of lateral dynamic displacement can be further mathematically modeled by harmonic superposition function.
(4)
By the comparison between simulative and monitoring lateral displacement effect from September 1 to September 3, mathematical models are feasible and effective to simulate time series of lateral girder displacement from main girder of Sutong Cable-Stayed Bridge.



---

*Source: 723152-2014-06-19.xml*