Non-Stationary Turbulent Wind Field Simulation of Long-Span Bridges Using the Updated Non-Negative Matrix Factorization-Based Spectral Representation Method

: Numerical simulation of the turbulent wind ﬁeld on long-span bridges is an important task in structural bu ﬀ eting analysis when it comes to the system non-linearity. As for non-stationary extreme wind events, some e ﬀ orts have been paid to update the classic spectral representation method (SRM) and the fast Fourier transform (FFT) has been introduced to improve the computational e ﬃ ciency. Here, the non-negative matrix factorization-based FFT-aided SRM has been updated to generate not only the horizontal non-stationary turbulent wind ﬁeld, but also the vertical one. Speciﬁcally, the evolutionary power spectral density (EPSD) is estimated to characterize the non-stationary feature of the ﬁeld-measured wind data during Typhoon Wipha at the Runyang Suspension Bridge (RSB) site. The coherence function considering the phase angles is utilized to generate the turbulent wind ﬁelds for towers. The simulation accuracy is validated by comparing the simulated and target auto- / cross-correlation functions. Results show that the updated method performs well in generating the non-stationary turbulent wind ﬁeld. The obtained wind ﬁelds will provide the research basis for analyzing the non-stationary bu ﬀ eting behavior of the RSB and other wind-sensitive structures in adjacent regions.


Introduction
Flexible long-span bridges are susceptible to extreme wind. As one of the wind-induced structural responses, bridge buffeting has been the concern of engineers for a long time [1]. The buffeting response, which is induced by turbulent winds, can be analyzed in both time [2] and frequency domain [3]. It is noted that the linear hypothesis is the prerequisite of the buffeting analysis in the frequency domain. Since the aerodynamic and geometrical effects may cause flexibilities and non-linear behaviors in long-span bridges, frequency-domain buffeting analysis is limited to most of long-span bridges [4][5][6]. As for time domain analysis, numerical integration is selected to solve the buffeting responses, hence, the non-linear and aerodynamic behaviors of long-span bridges can be considered [7,8]. To conduct time domain analysis, it is significant to accurately simulate the turbulent wind fields of wind-sensitive structures, since wind-induced responses are sensitive to the simulated wind loads [9].
The spectral representation method (SRM), which is an accurate and simplified method, is widely employed to generate the turbulent wind fields in structural wind engineering [5,10]. According to the stationary hypothesis, the turbulence is typically regarded as a zero-mean stationary stochastic process over a given time, which is determined by the power spectral density (PSD) [11,12]. However, field measurements indicate that extreme wind shows strong non-stationary features [13,14]. Therefore, The Runyang Yangtze River Bridge connects Yangzhou city and Zhenjiang, and is composed of the Runyang Suspension Bridge (RSB) (south side) and the Runyang Cable-stayed Bridge (north side). The main span of the RSB (see Figure 1) is 1490 m. Each two side spans are 470 m. The RSB is a hinged and simply supported steel-box-girder bridge. As a national key engineering project, the RSB is located on the Southeast coast of China. The meteorological survey indicates that this region belongs to subtropical monsoon climate. The cold and dry northern wind from the mainland is the main strong wind in winter, while the typhoon dominates in summer. Every year, the RSB is attacked by typhoons. The ultrasonic anemometer was installed on south tower (210 m in elevation) to collect wind data using a cylindrical coordinate system for further study of the typhoon characteristics at the RSB site. In September 2007, Typhoon Wipha passed through Jiangsu Province, the 3D strong wind data was collected by the anemometer. The sampling frequency was set as 20 Hz. The detailed information about the field measurement can be found in Wang et al. [20]. The sampling number in time was 40960 for better use of the FFT. Data collected from 20:00:00 to 20:34:08 on 19 September were utilized to estimate the EPSD. Daubechies wavelet (Db10) with a decomposition level of 11 was selected to separate the time-varying mean wind speed [21].
The stationarity of the turbulent wind flows was verified by the run test method with a significance level of 0.05 [22]. The field measured wind data and its time-varying mean are shown in Figure 2a,c). The turbulence is shown in Figure 2b,d).

Evolutionary Power Spectral Density (EPSD) Estimation of the Non-Stationary Turbulence
The wavelet transform, which can capture both the low-frequency and high-frequency information using different scale windows, was selected to calculate the EPSDs of the measured turbulence in Figure 2 [23]. In addition, the filtered harmonic wavelet (FHW), which is the updated version of the generalized harmonic wavelet, was selected to estimate the EPSD herein [24]. The EPSD for FHW scheme can be given by:

Evolutionary Power Spectral Density (EPSD) Estimation of the Non-Stationary Turbulence
The wavelet transform, which can capture both the low-frequency and high-frequency information using different scale windows, was selected to calculate the EPSDs of the measured turbulence in Figure 2 [23]. In addition, the filtered harmonic wavelet (FHW), which is the updated version of the generalized harmonic wavelet, was selected to estimate the EPSD herein [24]. The EPSD S ω i , t k for FHW scheme can be given by: where c F(m,n),k is the FHW transform coefficients. Equation (1) defines the local spectrum in the time and frequency regions.
The EPSDs of the measured turbulence in Figure 2b,d were estimated using the FHW with a total 2000 of frequency segmentations [25], as shown in Figure 3. The reliability of the estimated EPSDs was validated by comparing the averaging EPSDs with the power spectral density (PSDs) calculated by the Pwelch method [26]. Here, the averaging EPSD can be obtained by averaging the EPSD over time.
The EPSDs of the measured turbulence in Figure 2b,d were estimated using the FHW with a total 2000 of frequency segmentations [25], as shown in Figure 3. The reliability of the estimated EPSDs was validated by comparing the averaging EPSDs with the power spectral density (PSDs) calculated by the Pwelch method [26]. Here, the averaging EPSD can be obtained by averaging the EPSD over time.  The estimated EPSDs in along-wind and vertical directions are shown in Figure 3a,c, respectively. For clarity, the EPSDs with frequency over 2 Hz are omitted due to the relatively small values in Figure 3a,c. It is seen that the non-stationary characteristics of the measured turbulence are captured properly by the FHW scheme. Meanwhile, the averaging EPSDs agree quite well with the PSDs in Figure 3b,d, which indicates that the estimated EPSDs are accurate and consistent with the analytical theory [27].  The estimated EPSDs in along-wind and vertical directions are shown in Figure 3a,c, respectively. For clarity, the EPSDs with frequency over 2 Hz are omitted due to the relatively small values in Figure 3a,c. It is seen that the non-stationary characteristics of the measured turbulence are captured properly by the FHW scheme. Meanwhile, the averaging EPSDs agree quite well with the PSDs in Figure 3b,d, which indicates that the estimated EPSDs are accurate and consistent with the analytical theory [27].

The Updated Non-Negative Matrix Factorization (NMF)-Based Fast Fourier Transform (FFT)-Aided Spectral Representation Method (SRM)
3.1. Simulation of Non-Stationary Process Using the Classic SRM S(ω, t) is the EPSD of a zero-mean non-stationary process with components x 1 (t), x 2 (t), · · · , x n (t), which can be expressed as [15]: Appl. Sci. 2019, 9, 5506 5 of 17 S ij (ω, t) = S ii (ω, t)S jj (ω, t)γ ij (ω) (4) where ω represents the circular frequency; S ii (ω, t) = S ii (−ω, t) represents the auto-EPSD; S ij (ω, t) = S * ij (−ω, t) represents the cross-EPSD ('*' denotes the complex conjugate); γ ij (ω) are complex coherence functions between x i (t) and x j (t). The coherence matrix Γ(ω, t) can be represented by: [15] The cross-correlation function can be determined by [15]: The cross-correlation matrix for non-stationary processes is determined by t and t + τ (τ is the time lag). To simulate the one-dimensional, n variate (1D-nV) non-stationary stochastic process , where T represents the transpose, the S(ω, t) should be decomposed at each time instant t, into the following product [15]: Cholesky's method can finish this decomposition, and H(ω, t) can be represented by: where H ii (ω, t) are real and positive functions, while H ij (ω, t) are generally complex. The diagonal elements H ii (ω, t) satisfy the following relation: Furthermore, H ij (ω, t) can be written in polar form as: where Re and Im represent the real and imaginary parts of a complex number, respectively. Once the S(ω, t) is decomposed using Equations (8) and (9), the components of the x(t) can be simulated with the following expression [15]: where ω u represents the cutoff frequency; N represents the number of discretized frequencies; ∆ω is the frequency resolution; and Φ ml represent the independent random phase angles with uniform distribution from 0 to 2π. To avoid aliasing, the time resolution ∆t must satisfy the following relation: where M represents the discretized time number. It is pointed out that the simulation efficiency is unacceptable for engineering practice due to the time-dependent H(ω, t) in the double summation in Equation (13). Therefore, it is necessary to improve the simulation efficiency of the traditional non-stationary SRM.

The Updated Algorithm Considering Vertical Turbulent Wind Field
The essence of the non-negative matrix factorization (NMF) is a bound-constrained minimization problem. As for a non-negative matrix L ∈ R M×N and a pre-specified positive integer P < min{M, N}, NMF can be utilized to find two non-negative matrices W ∈ R M×P and V ∈ R P×N satisfying [28]: where WV is the NMF of L. The conventional approach proposed by Berry et al. [29] can be utilized to find the NMF of L: It is noted that Equation (19) is a standard optimization problem with bound constraint. The simple iteration algorithm proposed by Lee and Seung [30] is the widely-used approach to solve Equation (19).
At every iteration, W ip and V pj are multiplied by certain factors. When the iteration is stopped, W and V are strictly positive. It is noted that Equation (19) will become invariant if the stationary points of the W and V are reached in Euclidean distance [30].
The H jm (ω, t) can be further factorized using the NMF, then the FFT technique can be utilized in the non-stationary SRM.
Appl. Sci. 2019, 9, 5506 7 of 17 Substituting Equation (21) into Equation (13) and using Euler's formula satisfying: Note that the simulation efficiency of Equation (22) is much higher when the coherence matrix Γ(ω) is time-invariant. From Equation (3), S(ω, t) can be rewritten as [31]: It is obvious that Γ(ω) is a Hermitian matrix, which also has the positive property. When the elements in Γ(ω) are complex, the EPSD matrix can be rewritten in polar form [32]: Note that the phase angle in Equation (12) was the function coupled with time and frequency. However, it is only determined by frequency with regard to Equation (23). Thus, Θ can be expressed as: where ⊗ denotes the element-wise multiplication operator between two matrices. Similarly, Γ(ω) can be given by: where Γ(ω) represent the lagged coherence matrix [33]. According to Equation (8), S(ω, t) can be factorized as: where Λ(ω, t) is lower triangular matrix. The following relations have been proved when the phases admit the requirement that θ jm (ω) = θ jk (ω) − θ km (ω) [32]: Similar to Equation (23), S(ω, t) can be rewritten as: where B(ω) is the lower triangular matrix derived from Γ(ω) through Cholesky decomposition. That is: where According to Equations (30) and (33), the following relationship can be obtained: Based on Equations (31), (32) and (36), the components of the x(t) can be generated using: Taking advantage of the NMF, the square root of S ii (ω, t) can be factorized and given as: Substituting Equation (38) into Equation (37), the updated NMF-based non-stationary SRM can be rewritten as: Or It is clear that the FFT technique can be introduced for the second summation in Equation (40), which can enhance the simulation efficiency significantly. In addition, various auto-EPSDs and the complex coherence functions can be considered in Equation (40) for the simulation of both horizontal and vertical non-stationary turbulent wind fields for engineering structures.

The Hypotheses
Turbulent wind field simulation for long-span bridges is complicated due to the huge number of simulated points. According to the structural features of the suspension bridge, some hypotheses can be construed to enhance the simulation. In this work, non-stationary winds are regarded as a time-varying mean wind plus non-stationary turbulence: where u is the turbulence in the across-bridge direction, v is the turbulence in the along-bridge direction x and w is the turbulence in vertical direction Z, t is the time. U(Z, t) represents the time-varying mean wind. This is assuming that the EPSD matrices are uniform for all simulated points [34]. To simplify the calculation, the wind fields should be simulated in finite simulated points. Usually, only the correlation in spatial distributions of turbulent winds is considered in practice. Hence, a 3D-nV stochastic vector process can be simplified as three independent 1D-nV stochastic processes [35]. Specifically, turbulence v is not considered for the main girder since it is a linear component. As along-bridge vibrations of the towers are determined by the main cables, turbulence v for towers is not simulated. Meanwhile, turbulence w for towers is also not simulated [35]. In general, four independent 1D-nV stochastic processes can be assumed for the main girder and towers of the RSB based on the above hypotheses, listed in Table 1. Right tower u 30 7 The simulated points in the RSB is shown in Figure 4. There are 19 simulated points (from G1 to G19) uniformly distributed along the main girder (69.3 m in elevation), with the interval being 80.5 m. Similarly, there are 7 simulated points distributed uniformly along each tower. Specifically, T1 to T7 points in the left tower and T8 to T14 points in the right tower, with the interval being 30 m.

NMF of the Measured Decomposed EPSDs
The square root of the measured EPSDs in Figure 3a,c, can be factorized into the latent matrices W and V with the NMF. The integer P in Equation (38) is taken as 8, while the iteration tolerance is taken as 10 −8 in this factorization. The NMF results of the measured EPSDs are shown in Figure 5 with the target ones.

NMF of the Measured Decomposed EPSDs
The square root of the measured EPSDs in Figure 3a,c, can be factorized into the latent matrices W and V with the NMF. The integer P in Equation (38) is taken as 8, while the iteration tolerance is taken as 10 −8 in this factorization. The NMF results of the measured EPSDs are shown in Figure 5 with the target ones.
As shown in Figure 5, the NMF of the decomposed EPSD generally fit well with the targets in along-wind and vertical directions. Furthermore, the spectra snapshots of the NMF results at t = 200, 500, 1200 and 1500 s are compared with the targets. As shown in Figures 6 and 7, the spectrum snapshots agree well with the targets at the specific time instants, indicating that the NMF is an effective method to factorize the measured decomposed EPSDs.

NMF of the Measured Decomposed EPSDs
The square root of the measured EPSDs in Figure 3a,c, can be factorized into the latent matrices W and V with the NMF. The integer P in Equation (38) is taken as 8, while the iteration tolerance is taken as 10 −8 in this factorization. The NMF results of the measured EPSDs are shown in Figure 5 with the target ones. As shown in Figure 5, the NMF of the decomposed EPSD generally fit well with the targets in along-wind and vertical directions. Furthermore, the spectra snapshots of the NMF results at t = 200, 500, 1200 and 1500 s are compared with the targets. As shown in Figures 6 and 7, the spectrum snapshots agree well with the targets at the specific time instants, indicating that the NMF is an effective method to factorize the measured decomposed EPSDs.

Turbulent Wind Field Simulation for the Main Girder
The Davenport exponential function model is adopted to be the coherence function for the turbulence simulation of main girder [3]:

Turbulent Wind Field Simulation for the Main Girder
The Davenport exponential function model is adopted to be the coherence function for the turbulence simulation of main girder [3]: In Equations (42) and (43), λ represent the decay factor; ∆ is the distance between the successive simulated points; U(Z, t) represents the mean wind speed, which can be obtained by averaging the U(Z, t). The mean wind speed profile at the RSB site is consistent with the power law: where U 1 and U 2 denote wind speeds at altitude Z 1 and Z 2 , respectively, and α is a dimensionless exponent which depends on roughness of terrain, taken as measured value 0.12 herein [35]. Based on Equation (36), the explicit form of the lagged decomposed EPSD matrix Λ can be expressed as [36]: Here, the Davenport coherence function is adopted and the decay factor is taken as seven. Table 2 lists the detailed information in the simulation. As shown in Figure 8, 5000 samples (both in Along-wind and Vertical direction) are simulated using the measured EPSDs. Then, the auto-/cross-correlation functions are estimated using the generated samples. The generated wind fields are validated by comparing the estimated correlation functions and the target ones at various time differences τ, as shown in Figure 9.
Here, the Davenport coherence function is adopted and the decay factor is taken as seven. Table  2 lists the detailed information in the simulation. As shown in Figure 8, 5000 samples (both in Along-wind and Vertical direction) are simulated using the measured EPSDs. Then, the auto-/cross-correlation functions are estimated using the generated samples. The generated wind fields are validated by comparing the estimated correlation functions and the target ones at various time differences τ, as shown in Figure 9.      Figure 8 shows that the generated wind samples in along-wind and vertical directions exhibit prescribed non-stationary features. In Figure 9, the estimated auto-correlation functions fit well with the target one, while the discrepancy between the cross-correlation function and the target one is relatively wide. This phenomenon will be improved if many more samples are considered.

Turbulent Wind Field Simulation for Towers
As for towers, the mean wind speed approximately obeys the distribution of the exponential relationship in Equation (44). The following empirical coherence function is selected to generate the non-stationary turbulence for towers in the along-wind direction [32].
According to Equation (46), phase angles θ jm (ω) in Equation (32) can be determined by: In the wind field, the apparent wave velocity v app can be modified as [37]: where U j and U m denote the mean wind speed at different elevation; C θ represents the experimental appropriate coefficient [3]. Obviously, v app is varying for the simulation points of the towers with different mean wind speeds. However, it can be approximated to be a constant by calculating the mean wind speeds on two end points, and for that the variation is usually acceptable [32]. Table 3 lists the detailed information in the simulation. As shown in Figure 10, the non-stationary turbulent wind fields of towers in the along-wind direction are shown. For comparison, 5000 samples are utilized to estimate the ensemble auto-/cross-correlation functions at various time differences τ. The estimated auto-/cross-correlation functions and corresponding target ones are shown in Figure 11. In the wind field, the apparent wave velocity app v can be modified as [37]: where j U and m U denote the mean wind speed at different elevation; C θ represents the experimental appropriate coefficient [3]. Obviously, app v is varying for the simulation points of the towers with different mean wind speeds. However, it can be approximated to be a constant by calculating the mean wind speeds on two end points, and for that the variation is usually acceptable [32]. Table 3 lists the detailed information in the simulation.  Figure 10, the non-stationary turbulent wind fields of towers in the along-wind direction are shown. For comparison, 5000 samples are utilized to estimate the ensemble auto-/crosscorrelation functions at various time differences τ. The estimated auto-/cross-correlation functions and corresponding target ones are shown in Figure 11.    As shown in Figure 10, the non-stationary turbulent wind fields of towers are constructed from the measured EPSD in along-wind direction with the coherence function in Equation (46). It is seen that the estimated correlation functions generally fit well with the target ones in Figure 11, indicating the reliability of the generated non-stationary turbulent wind fields of the towers. As shown in Figure 10, the non-stationary turbulent wind fields of towers are constructed from the measured EPSD in along-wind direction with the coherence function in Equation (46). It is seen that the estimated correlation functions generally fit well with the target ones in Figure 11, indicating the reliability of the generated non-stationary turbulent wind fields of the towers.

Conclusions
The existing FFT-aided NMF-based SRM, which cannot simulate the vertical turbulent wind field, has been updated to simulate a non-stationary turbulent wind field in both horizontal and vertical directions. NMF can be utilized to decouple the EPSD matrix, the FFT technique can be then utilized to improve the simulation efficiency of the non-stationary SRM. The turbulent wind fields of the RSB are generated based on the field-measured EPSDs during typhoon Wipha. The non-stationary turbulent wind fields for towers are generated considering the coherence function with the phase angles. Relative comparisons are conducted to validate the reliability and accuracy of the generated wind fields. The following conclusions can be drawn accordingly: • The existed NMF-based SRM for horizontal non-stationary turbulent wind field simulation is updated to also generate the vertical non-stationary turbulent wind field. With the aid of the FFT, the simulation efficiency is enhanced.

•
In the NMF-based SRM, the phase angles can be considered in the coherence functions. The complex EPSD matrix can be transformed into the real modulus matrix after the phase was separated. With this treatment, the efficiency of decomposition can be improved.

•
The estimated auto-/cross-correlation functions generally fit well with the target ones, validating the effectiveness and reliability of the non-stationary turbulent wind fields generated. The presented wind field simulation method can provide research basis for non-stationary buffeting analysis of the RSB.