# Evolutionary Spectral Analyses of a Powerful Typhoon at the Sutong Bridge Site Based on the HHT

Authors: Lin Ma; Da-jun Zhou; Ai-min Yuan; Hao Wang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

To investigate the nonstationary characteristics of strong typhoons, this paper considers the evolutionary spectral characteristics of strong typhoons based on the Hilbert-Huang transform (HHT). Discrete expressions are determined for the evolutionary spectral analysis based on the HHT. The study indicates that the classic empirical mode decomposition (EMD) method fails to extract all of the high-frequency fluctuations from wind velocity data, and the time-averaged power spectrum obtained directly using the HHT cannot provide the true wind velocity spectrum. The degrees of nonstationarity of different-order IMF components are analysed and a synthesized method of analysing the evolutionary spectrum and time-averaged power spectrum of a strong typhoon is proposed. To avoid the energy leakage problem that exists in HHT spectral analyses, the Gram-Schmidt method is used to orthogonalize the intrinsic mode function (IMF) components. The study indicates that when the orthogonalization is implemented in accordance with the sequence from high-order IMF components to low-order ones, the orthogonalized components retain the same good Hilbert property as that of the IMFs. The synthesized method proposed yields a time-averaged power spectrum that is consistent with the Fourier spectrum in value and can produce the energy distributions of a typhoon in the time and frequency domains simultaneously.


---

## Body

## 1. Introduction

As the span of a bridge increases, the wind load plays an increasingly important role in its design. Determining an accurate wind load as a design parameter is critical for ensuring the safety of an engineering structure. The major wind characteristics at the site of a bridge or building include the average wind speed, wind direction, wind profile, turbulence intensity, gust factor, turbulence integral scale, and power spectral density. The power spectral density is an important parameter of fluctuating wind and is used to describe the energy distribution of the turbulence over different frequencies. Currently, the fluctuating wind spectra that are most widely used in wind resistance design for bridges or buildings include the Davenport [1], Harris [2], Kaimal et al. [3], and Teunissen [4] spectra, among others.The literature on wind turbulence characteristics is extensive. Lots of authors have published the results of full-scale measurements and wind tunnel tests, whereas others have proposed new mathematical methods of describing turbulence characteristics or algorithms for the numerical simulation of wind fields. On the topic of theoretical models, Solari [5] published a critical review of the models available for representing the longitudinal component of atmospheric turbulence. Starting from a critical review of the state of the art, Solari and Piccardo [6] proposed a unified model of atmospheric turbulence that is especially well suited to determining the 3D gust-excited response of a structure. Unlike in classic models, all parameters addressed in [6] are assigned based on first- and second-order statistical moments derived from a broad set of selected experimental measurements. Solari and Tubino [7] investigated the two-point coherence function of the different components of turbulence and suggested a physical principle to establish an appropriate model of the two-point coherence function of the longitudinal and vertical components, thereby completing the statistical model of turbulence. The proper orthogonal decomposition (POD) approach offers efficient tools for formulating a model of a turbulence field based on principal components. Tubino and Solari [8] proposed a representation of turbulence based on a double application of the POD.In recent years, researchers have performed many field measurements of the wind characteristics at bridge sites [9–15]. Using the wind data obtained from the Structural Health Monitoring System (SHMS), Xu et al. [9] studied the mean and turbulent wind characteristics of Typhoon Victor at the site of the Tsing Ma Bridge in Hong Kong. It was found that certain features of the wind structure and bridge response were difficult to consider in the analytical process that is currently used to predict the buffeting response of long suspension bridges, as the bridge is surrounded by a complex topography and the wind direction of Typhoon Victor varied during its passage. Miyata et al. [10] analysed full-scale measurement data observed on the Akashi-Kaikyo Bridge during strong typhoons. The power spectral density and spatial correlation of the longitudinal velocity fluctuations were analysed. These authors found that the spatial correlation of two points was well represented not by an exponential formula but rather by an alternative coherence function based on isotropic turbulence theory. To derive a model of turbulence suitable for buffeting calculations regarding the Stonecutters Bridge in Hong Kong, Hui et al. [11, 12] studied the mean wind, turbulence intensities, wind power spectra, integral length scales, and wind coherences of the bridge.Using statistical theory, Wang et al. [13] analysed the strong wind characteristics of Typhoon Matsa, the buffeting response characteristics of the cable and deck of the Runyang Yangtze Bridge, and the variation in the buffeting response RMS versus the wind speed. The results obtained in this study served to validate the credibility of current techniques in buffeting response analysis theory. Based on field measurements of Typhoon Nuri at the Macao Friendship Bridge site, Song et al. [14] proposed that the value of the integral scale increased when the eye wall of Typhoon Nuri passed over the field measurement site. In the eye wall of the typhoon, the horizontal spatial correlation was relatively strong, and the horizontal spatial correlation spectrum decayed slowly with increasing frequency. Liu et al. [15] compared the wind characteristics of typhoons and strong monsoons at the site of the Xi-Hou-Men Bridge. The similarities and differences of wind characteristics between typhoons and monsoons were analysed. Wang et al. [16] analysed the recorded real-time wind data measured at the Sutong Yangtze Bridge site in detail to generate the wind-rose diagram, mean wind speed and direction, turbulence intensity, turbulence integral scale and power spectral density, and conducted comparative analyses among the inhomogeneous wind characteristics of three strong wind events, including the Northern wind, Typhoon Kalmaegi and Typhoon Fung-Wong. Wang et al. [17] analysed the buffeting responses of the cable-stayed Sutong Bridge with the wind spectra used in the design phase and those obtained experimentally from a long-term SHM system to compare the actual buffeting response with design predictions.In previous studies, the Fourier transform has often been used to obtain wind spectrum information based on the assumption of stationarity. However, several studies have shown that powerful typhoons exhibit strongly nonstationary characteristics and that lots of abrupt, instantaneous pulses occur in strong typhoons. Therefore, to accurately describe the nonstationary characteristics of powerful typhoons, it is necessary to study their evolutionary power spectra [18].Because of the limitations of theory, it is difficult to use the Fourier transform to analyse a nonstationary signal. First, Fourier theory defines harmonic components over all time, and each frequency exists in a complete harmonic form in the time domain; thus, many additional harmonic components are required to simulate uneven nonstationary data, which may induce false harmonic components and energy divergence. Second, the Fourier transform lacks the capability of simultaneous positioning in both the time and frequency domains. The spectral analysis results obtained using the Fourier transform represent the time-averaged spectral distribution. Finally, the time resolution and frequency resolution of a signal are related to each other, and the Fourier transform fails to automatically adjust the resolutions in the time and frequency domains depending on the signal characteristics.To address the inadequacy of Fourier analysis, researchers have proposed time-frequency analysis methods in which a time-frequency transform is applied to expand the signal energy in the time-frequency plane to reflect the process of variation of a signal in the time and frequency domains. Such a time-frequency analysis can describe the changing regularity of a signal’s spectrum over time. Commonly used time-frequency analysis methods include the short-time Fourier transform (STFT), the Winger-Ville distribution, the wavelet transform, and the Hilbert-Huang transform (HHT).Huang et al. [19] proposed a revolutionary method of signal decomposition called empirical mode decomposition (EMD). The EMD components possess the Hilbert property, that is, a time-frequency distribution of these components that bears an actual physical meaning can be obtained by applying the Hilbert transform. The US National Aeronautics and Space Administration refers to this new time-frequency analysis method as the HHT. The essence of the method is to decompose an original signal into its different fluctuations and trends at different scales, thereby obtaining a series of data with different size characteristics.Although researchers have made considerable progress in the field of time-frequency analysis, the amplitude spectrum obtained via the direct wavelet transform or the HHT has a completely different mathematical definition from and offers no quantitative comparability with the traditional Fourier spectrum. It is difficult to explain the physical meaning of the amplitude spectrum obtained via a wavelet transform or the HHT. Priestley [20] proposed the concept of an evolutionary spectrum by generalizing the power spectrum of a stationary random process. The evolutionary spectrum, as a conceptual extension of the power spectrum of a stationary random process, can accurately describe the changing characteristics of a spectrum over time. The evolutionary spectrum at any instant of time has the same physical significance as the power spectrum of a stationary random process. Therefore, the evolutionary spectrum is considered to possess a good physical interpretation. However, parameter estimation for an evolutionary spectrum has always been a difficult problem, in which both the power spectral density function and the time and frequency modulation functions must be estimated.Spanos and Failla [21] studied parameter estimation for an evolutionary spectrum based on the wavelet transform. They calculated the evolutionary spectrum using the frequency-domain expression of a wavelet function. Ding et al. [22] discussed the estimation of the evolutionary spectrum of a nonstationary process using the time-domain wavelet expression and wavelet transform coefficient and studied the evolutionary spectra of the wind velocity and bridge response when Typhoon Matsa crossed the Runyang Yangtze Bridge site.Wen and Gu [18] considered the HHT spectral analysis of a nonstationary process and presented the conversion relationship between the Hilbert marginal spectrum and the Fourier energy spectrum. However, that work did not consider the energy leakage problem caused by incomplete orthogonality of the EMD components. Hu and Chen [23] proposed a method of estimating the local spectral density based on orthogonal EMD components, applied the method to the acceleration records of the El Centro and Taft waves, and studied the evolutionary spectral characteristics of these seismic waves.In the literature [18, 23], it is assumed that the energy of each component is always equal to half the energy of the corresponding analytical signal. As noted by Wen and Gu [18], for a high-frequency component, this assumption leads to only a small error because a sufficient number of waves are present. However, this assumption is not always valid for a low-frequency component. When the amplitude of a low-frequency component is large, a large error could be incurred. Furthermore, the above-mentioned work failed to discuss the problem of mode mixing in EMD, which leads to a failure of HHT analyses to provide complete power spectrum information. The Hilbert energy spectrum presented in [18] is consistent with the Fourier power spectrum at low frequencies, but it is considerably lower than the Fourier power spectrum at high frequencies. However, regardless of the analysis method used, the total energy spectrum or time-averaged power spectrum should be consistent or should at least satisfy the conservation of energy. Therefore, it appears that there is still some deficiency in the current method of HHT evolutionary spectrum analysis.This paper considers the evolutionary spectrum characteristics of a powerful typhoon based on the HHT. Discrete expressions for the evolutionary spectrum based on the HHT are deduced, and the Gram-Schmidt method is used to orthogonalize the EMD components to avoid energy leakage during the analysis of the typhoon’s evolutionary spectrum. The paper discusses the details of the problem of the missing high-frequency fluctuations in the time-averaged spectrum obtained using the HHT. Based on analyses of the degrees of nonstationarity of the EMD components, this paper proposes a synthesized method for studying the evolutionary spectrum and time-averaged spectrum of a typhoon. This synthesized method can avoid the deficiency of the HHT when analysing the evolutionary spectrum of a strong typhoon and can be used to study the nonstationarity characteristics of powerful typhoons.


## 2. Brief Review of the Hilbert-Huang Transform

### 2.1. Instantaneous Frequency

The widely accepted definition of the instantaneous frequency is that it is the derivative of the phase of the signal [24–26]. To obtain the instantaneous frequency, a real signal x(t) should be converted into an analytic signal z(t) in the form of a complex function. The corresponding imaginary part can be obtained from the Hilbert transform of the original signal: (1)yt=1πP∫-∞∞xτt-τdτ,where P denotes the Cauchy principal value.From this definition, the following analytic signal is obtained:(2)zt=xt+iyt.Although there are many ways to define the imaginary part of the signal, the Hilbert transform provides a means of defining the imaginary part such that the result is unique. Equation (2) can be expressed in polar form as(3)zt=ateiθt,where(4)at=x2t+y2t1/2,θt=arctan⁡ytxt.In this case, the instantaneous frequency can be defined as(5)ωt=dθtdt.The polar form of the analytic signal reflects the physical meaning of the Hilbert transform. It represents the time-varying amplitude and frequency of a signal. However, the Hilbert transform is only applicable to a monocomponent signal, namely, either a signal with a single frequency component or a narrow-band signal. For a multicomponent signal, the instantaneous frequency defined by the Hilbert transform is not necessarily a single-valued function of time [26]. Because of the lack of a precise definition of “monocomponent,” the “narrow-band” requirement has been adopted as the limitation placed on the data to ensure that the instantaneous frequency is meaningful. However, as noted by Huang et al. [19], this bandwidth limitation as defined in the global sense is overly restrictive yet it simultaneously lacks precision.


### 2.2. Empirical Mode Decomposition

Huang et al. [19] proposed the concept of the intrinsic mode function (IMF) to ensure the effectiveness of the Hilbert transform. An IMF is a function that satisfies two conditions: (1) in the entire dataset, the number of extrema and the number of zero crossings must either be equal or differ by at most one, and (2) at any point, the mean value of the envelope defined by the local maximum and that defined by the local minimum is zero. The first condition is similar to the traditional narrow-band requirements for a stationary Gaussian process. The second condition ensures that the instantaneous frequency will not exhibit unwanted fluctuations induced by asymmetric wave forms. Huang et al. [19] noted that even under the worst conditions, the instantaneous frequency defined by the Hilbert transform of an IMF is still consistent with the physics of the system under study. The IMF concept offers a more effective method of judging whether a signal is monocomponent. In addition, an IMF is not restricted to being a narrow-band stationary signal. An IMF can be modulated in both amplitude and frequency and can be a nonstationary signal.Generally, most real signals are not IMFs. A real signal must be decomposed into IMF components. Huang et al. [19] innovatively proposed the EMD method to address both nonstationary and nonlinear data. A specific algorithm for EMD can be found in the relevant literature. The result of the decomposition takes the following well-known form:(6)xt=∑j=1nCjt+rnt,where Cj(t) (j=1,2,…,n) denotes the jth component of the original signal and rn(t) is the residual signal. An IMF represents an intrinsic fluctuation of a signal, and each mode function at a different time scale corresponds to a wave pattern in the original signal. Because IMFs are extracted from a signal itself, they have self-adaptive features.A signal should be sifted many times before a qualified IMF component can be acquired. Unfortunately, infinitely repeated sifting could obliterate the physically meaningful amplitude modulations, making the resulting IMF a pure frequency-modulated signal of constant amplitude. To guarantee that the IMF components retain both amplitude and frequency modulations, it is necessary to determine a suitable termination criterion for the sifting process. To accomplish this task, Huang et al. [19] proposed a method of limiting the size of the standard deviation of two consecutive sifting results. They proposed that the reasonable standard deviation threshold lies in the range of 0.2~0.3. As an improvement to this criterion, Rilling et al. [27] proposed a criterion basedon 2 thresholds, θ1 and θ2, with the intent of guaranteeing globally small fluctuations while accounting for locally large excursions. They defined the mode amplitude as a(t)=(emax(t)-emin(t))/2, the mean as m(t)=(emax(t)+emin(t))/2, and the evaluation function as σ(t)=|m(t)/a(t)|. The functions emax(t) and emin(t) denote the upper and lower envelope lines, respectively. The sifting is iterated until σ(t)<θ1 for some prescribed fraction (1-α) of the total duration, whereas the criterion is σ(t)<θ2 for the remaining fraction. The threshold values suggested by Rilling et al. are as follows: α=0.05, θ1=0.05 and θ2=10θ1.The ending point effect is an important problem in the EMD algorithm. Because the two endpoints of the signal or the data are not always the extrema, the data at the end points may not be completely encompassed by the upper and lower envelopes. The error induced by the endpoint effect will spread inward during repeated sifting processes and “contaminate” the entire data series; under severe conditions, the endpoint effect will make the components meaningless. In this paper, the mirror method proposed by Rilling et al. [27] is adopted to mitigate the endpoint effect. In this method, extremum points are added by mirroring the extrema adjacent to the end points. The explicit steps of this method can be found in the literature [28].


### 2.3. Hilbert-Huang Transform

By applying the Hilbert transform to each IMF component, one can obtain the transient spatial spectrum of the components. The transient spatial spectrum of the original signal can be obtained by superposing of all those components:(7)Xt=Re∑jnajtei∫ωjtdt,where Re denotes the real part of a complex number. In (7), X(t) is not only a function of time but also a function of frequency. In this way, the signal is expanded over the time-frequency plane. The time-frequency distribution of the amplitude aj(t) is designated as the Hilbert amplitude spectrum, H(ω,t), or simply the Hilbert spectrum. This method of transient spatial spectrum analysis is called the Hilbert-Huang transform (HHT). The HHT spectrum describes the regularity of variation in a signal’s transient amplitude and frequency over time. In the Hilbert spectrum, a transient frequency represents a wave of such a frequency that occurs locally in time. This method is different from Fourier analysis, in which each frequency exists in a complete harmonic form over all time.


## 2.1. Instantaneous Frequency

The widely accepted definition of the instantaneous frequency is that it is the derivative of the phase of the signal [24–26]. To obtain the instantaneous frequency, a real signal x(t) should be converted into an analytic signal z(t) in the form of a complex function. The corresponding imaginary part can be obtained from the Hilbert transform of the original signal: (1)yt=1πP∫-∞∞xτt-τdτ,where P denotes the Cauchy principal value.From this definition, the following analytic signal is obtained:(2)zt=xt+iyt.Although there are many ways to define the imaginary part of the signal, the Hilbert transform provides a means of defining the imaginary part such that the result is unique. Equation (2) can be expressed in polar form as(3)zt=ateiθt,where(4)at=x2t+y2t1/2,θt=arctan⁡ytxt.In this case, the instantaneous frequency can be defined as(5)ωt=dθtdt.The polar form of the analytic signal reflects the physical meaning of the Hilbert transform. It represents the time-varying amplitude and frequency of a signal. However, the Hilbert transform is only applicable to a monocomponent signal, namely, either a signal with a single frequency component or a narrow-band signal. For a multicomponent signal, the instantaneous frequency defined by the Hilbert transform is not necessarily a single-valued function of time [26]. Because of the lack of a precise definition of “monocomponent,” the “narrow-band” requirement has been adopted as the limitation placed on the data to ensure that the instantaneous frequency is meaningful. However, as noted by Huang et al. [19], this bandwidth limitation as defined in the global sense is overly restrictive yet it simultaneously lacks precision.


## 2.2. Empirical Mode Decomposition

Huang et al. [19] proposed the concept of the intrinsic mode function (IMF) to ensure the effectiveness of the Hilbert transform. An IMF is a function that satisfies two conditions: (1) in the entire dataset, the number of extrema and the number of zero crossings must either be equal or differ by at most one, and (2) at any point, the mean value of the envelope defined by the local maximum and that defined by the local minimum is zero. The first condition is similar to the traditional narrow-band requirements for a stationary Gaussian process. The second condition ensures that the instantaneous frequency will not exhibit unwanted fluctuations induced by asymmetric wave forms. Huang et al. [19] noted that even under the worst conditions, the instantaneous frequency defined by the Hilbert transform of an IMF is still consistent with the physics of the system under study. The IMF concept offers a more effective method of judging whether a signal is monocomponent. In addition, an IMF is not restricted to being a narrow-band stationary signal. An IMF can be modulated in both amplitude and frequency and can be a nonstationary signal.Generally, most real signals are not IMFs. A real signal must be decomposed into IMF components. Huang et al. [19] innovatively proposed the EMD method to address both nonstationary and nonlinear data. A specific algorithm for EMD can be found in the relevant literature. The result of the decomposition takes the following well-known form:(6)xt=∑j=1nCjt+rnt,where Cj(t) (j=1,2,…,n) denotes the jth component of the original signal and rn(t) is the residual signal. An IMF represents an intrinsic fluctuation of a signal, and each mode function at a different time scale corresponds to a wave pattern in the original signal. Because IMFs are extracted from a signal itself, they have self-adaptive features.A signal should be sifted many times before a qualified IMF component can be acquired. Unfortunately, infinitely repeated sifting could obliterate the physically meaningful amplitude modulations, making the resulting IMF a pure frequency-modulated signal of constant amplitude. To guarantee that the IMF components retain both amplitude and frequency modulations, it is necessary to determine a suitable termination criterion for the sifting process. To accomplish this task, Huang et al. [19] proposed a method of limiting the size of the standard deviation of two consecutive sifting results. They proposed that the reasonable standard deviation threshold lies in the range of 0.2~0.3. As an improvement to this criterion, Rilling et al. [27] proposed a criterion basedon 2 thresholds, θ1 and θ2, with the intent of guaranteeing globally small fluctuations while accounting for locally large excursions. They defined the mode amplitude as a(t)=(emax(t)-emin(t))/2, the mean as m(t)=(emax(t)+emin(t))/2, and the evaluation function as σ(t)=|m(t)/a(t)|. The functions emax(t) and emin(t) denote the upper and lower envelope lines, respectively. The sifting is iterated until σ(t)<θ1 for some prescribed fraction (1-α) of the total duration, whereas the criterion is σ(t)<θ2 for the remaining fraction. The threshold values suggested by Rilling et al. are as follows: α=0.05, θ1=0.05 and θ2=10θ1.The ending point effect is an important problem in the EMD algorithm. Because the two endpoints of the signal or the data are not always the extrema, the data at the end points may not be completely encompassed by the upper and lower envelopes. The error induced by the endpoint effect will spread inward during repeated sifting processes and “contaminate” the entire data series; under severe conditions, the endpoint effect will make the components meaningless. In this paper, the mirror method proposed by Rilling et al. [27] is adopted to mitigate the endpoint effect. In this method, extremum points are added by mirroring the extrema adjacent to the end points. The explicit steps of this method can be found in the literature [28].


## 2.3. Hilbert-Huang Transform

By applying the Hilbert transform to each IMF component, one can obtain the transient spatial spectrum of the components. The transient spatial spectrum of the original signal can be obtained by superposing of all those components:(7)Xt=Re∑jnajtei∫ωjtdt,where Re denotes the real part of a complex number. In (7), X(t) is not only a function of time but also a function of frequency. In this way, the signal is expanded over the time-frequency plane. The time-frequency distribution of the amplitude aj(t) is designated as the Hilbert amplitude spectrum, H(ω,t), or simply the Hilbert spectrum. This method of transient spatial spectrum analysis is called the Hilbert-Huang transform (HHT). The HHT spectrum describes the regularity of variation in a signal’s transient amplitude and frequency over time. In the Hilbert spectrum, a transient frequency represents a wave of such a frequency that occurs locally in time. This method is different from Fourier analysis, in which each frequency exists in a complete harmonic form over all time.


## 3. Orthogonal EMD Components of the Wind Velocity Data of a Strong Typhoon

It is necessary to ensure that no energy leakage occurs in the evolutionary spectrum when the squares of the EMD components are used to calculate the energy of the original data. Thus, the IMF components should be orthogonal and complete. The components obtained through the EMD method proposed by Huang et al. [19] are approximately orthogonal. However, numerical examples have shown that the degree of orthogonality of the IMFs is not ideal. The orthogonality index, which includes the overall orthogonality index and the orthogonality index for any two components, is used to judge the degree of orthogonality:(8)IOT=∑j=1n+1∑k=1,k≠jn+1∫0Tcjtcktdt∫0TX2tdt=∑j=1n+1∑k=1,k≠jn+1∑i=1Ncjicki∑i=1NXi2,(9)IOjk=∫0Tcjtcktdt∫0Tcj2tdt+∫0Tck2tdt=∑i=1Ncjicki∑i=1Ncji2+cki2,where IOT and IOjk denote the overall index of orthogonality and the orthogonality index for any two components, respectively. When all IMFs are completely orthogonal, these indices are zero. In this case, the sum of the energies of all IMFs is equal to the total energy of the original data; that is, there is no energy leakage.By applying the Gram-Schmidt orthogonalization method to these IMFs, one can obtain completely orthogonal components. The explicit procedure is as follows.(1) By decomposing a signal using the EMD method, the IMFs ofX(t) are obtained, which are denoted by c1t,c2t,…,cnt. In addition, c-1t,c-2t,…,c-nt are used to denote the orthogonalized results of the IMF components. The first orthogonal component is obtained by setting c-1t=c1t.(2) As seen from the sifting process ofc2t, one cannot ensure that c2t is completely orthogonal to c-1t. To obtain the second orthogonal component of X(t), one should subtract the component relevant to c-1t from c2t:(10)c-2t=c2t-β21c-1t,where c-2t is the 2nd orthogonal component of X(t) and β21 is the coefficient of orthogonalization between c2t and c-1t. To obtain β21, one should multiply (10) by c-1t and integrate it over time. Considering the orthogonality of c-2t and c-1t, one can deduce coefficient β2 as follows: (11)∫0Tc-2tc-1tdt=∫0Tc2tc-1tdt-β21∫0Tc-12tdt=0,(12)β21=∫0Tc2tc-1tdt∫0Tc-12tdt.Below, (12) is expressed in discrete form: (13)β21=c2c-1Tc-1c-1T=∑i=1Nc2ic-1i∑i=1Nc-1i2.(3) Using the same approach, the(j+1)th orthogonal component can be expressed as follows:(14)c-j+1t=cj+1t-∑i=1jβj+1,ic-it,where c-j+1(t) (j=1,2,…,n) denotes the (j+1)th orthogonal component. By multiplying (14) by c-k(t) (k≤j) and integrating the equation over time, one obtains(15)∫0Tc-j+1tc-ktdt=∫0Tcj+1tc-ktdt-∑i=1jβj+1,i∫0Tc-itc-ktdt=0.From the orthogonality ofc-k(t) with respect to c-i(t) (i=1,2,…,j+1 and i≠k), one can deduce the coefficient βj+1,k as follows:(16)βj+1,k=∫0Tcj+1tc-ktdt∫0Tc-k2tdt.The discrete form of (16) is (17)βj+1,k=cj+1c-kTc-kc-kT=∑m=1Ncj+1,mc-k,m∑m=1Nc-k,m2.After the above procedures are applied, the original signal is expressed in the following form:(18)xt=∑j=1ncjt+rnt=∑j=1n∑i=1jβj,ic-it+rnt=∑i=1naic-it+rnt=∑i=1nci∗t+rnt,where(19)ci∗t=aic-iti=1,2,…,n,ai=∑j=inβj,ii=1,2,…,n,βj,i=1j=i.Because the componentsc-i(t) (i=1,2,…,n) satisfy the condition of complete orthogonality, the ci∗(t) (i=1,2,…,n) are also completely orthogonal to each other because a linear transform does not change the orthogonality relationship between components. In this way, the signal x(t) is decomposed into the sum of the orthogonal IMFs (OIMFs) ci∗(t) (i=1,2,…,n) and the residue rn(t).In this paper, the wind velocity of Typhoon Haikui at the Sutong Yangtze Bridge site is selected as the object of analysis. The sampling frequency of the measured wind data is 1 Hz. Figure1 shows the wind velocity data recorded at the top of the bridge’s south tower during an hour when the mean wind reached its maximum. The time history of the wind velocity data contains 3,544 data points; this value is less than 3,600 because a small amount of data appears to have been omitted in the data collection. The wind fluctuations are obtained by subtracting the 10 min mean wind from the original wind velocity data. The IMFs and the residual term obtained by applying the EMD method to the time history of the wind fluctuations are shown in Figure 2. The OIMF components obtained by orthogonalizing these components using the Gram-Schmidt method are presented in Figure 3.Figure 1
Wind velocity of “Haikui” at the top of the bridge’s south tower.
(a)
Original wind
(b)
Fluctuating windFigure 2
IMF components of the wind velocity of “Haikui.”Figure 3
Orthogonalized components of the wind velocity of “Haikui.”Figure2 shows that the wind fluctuations are decomposed into 9 IMFs and one residue. The instantaneous frequency of the first IMF component is the highest, and those of the other components sequentially decrease. From an energy perspective, the sum of the squares of the fluctuating wind data is 29581, and the sum of the squares of all OIMFs and the corresponding residue is the same; however, the sum of the squares of IMFs and the residue is 31274. This extra energy originates from the incomplete orthogonality of the IMFs. The orthogonality indices calculated using (8) are presented in Table 1. The indices calculated between the IMFs are shown in the upper triangular portion of the table, whereas those calculated between the OIMFs are given in the lower triangular portion. Table 1 shows that the orthogonal indices between the OIMFs are less than 10−17 and approximately equal to zero, which are far less than those between the IMFs. The overall index of orthogonality of the IMFs is −0.0572, whereas that of the OIMFs is nearly zero. Thus, it can be concluded that the Gram-Schmidt orthogonalization method can ensure the complete orthogonality of the components. When the sum of the squares of the OIMFs is used to calculate the energy, no energy leakage will be produced; this is the premise of the evolutionary spectrum analysis presented in subsequent sections.Table 1
Orthogonality index of the components.
(a)
IMF
c
1
c
2
c
3
c
4
c
5
c
1
0.5000
0.0195
-
0.0066
-
0.0081
-
0.0032
c
2
8
.
4
e
-
018
0.5
0.0161
-
0.0055
0.0026
c
3
-
1.2
e
-
018
2
.
6
e
-
018
0.5000
0.0309
-
0.0027
c
4
3
.
5
e
-
018
-
5
.
6
e
-
018
1
.
3
e
-
017
0.5000
-
0.0121
c
5
7.2
e
-
018
-
1.7
e
-
017
-
1.7
e
-
018
-
1.9
e
-
017
0.5000
c
6
-
5
.
1
e
-
018
6.3
e
-
018
-
2.9
e
-
019
2.9
e
-
018
1
.
7
e
-
018
c
7
-
9
.
3
1
e
-
018
-
4.6
e
-
018
2.5
e
-
018
-
4.7
e
-
018
-
3
.
8
e
-
018
c
8
8.3
e
-
018
7.1
e
-
018
1.1
e
-
018
-
8.0
e
-
019
-
5.6
e
-
019
c
9
-
5
.
6
e
-
018
1
.
5
e
-
018
-
1.4
e
-
018
4
.
1
e
-
019
4.4
e
-
018
r
9
8.7
e
-
019
-
1
.
2
e
-
018
1
.
1
e
-
019
-
4.7
e
-
019
-
6
.
7
e
-
019
(b)
IMF
c
6
c
7
c
8
c
9
r
9
c
1
0.0004
-
0.0008
0.0079
-
0.0068
0.0000
c
2
-
0.0059
-
0.0008
-
0.0023
0.0018
0.0049
c
3
0.0085
-
0.0011
-
0.0005
0.0038
-
0.0004
c
4
-
0.0040
0.0110
-
0.0030
0.0110
-
0.0028
c
5
0.0187
-
0.0244
0.0192
-
0.0150
-
0.0015
c
6
0.5000
-
0.0009
-
0.0110
-
0.0092
0.0002
c
7
5
.
2
e
-
018
0.5000
-
0.0252
0.0183
-
0.0127
c
8
-
1
.
5
e
-
018
-
5
.
4
e
-
018
0.5000
-
0.1445
0.0059
c
9
-
4
.
3
e
-
018
1.3
e
-
017
6.0
e
-
017
0.5000
-
0.0342
r
9
-
2.2
e
-
018
6
.
1
e
-
018
-
3.1
e
-
018
8
.
5
e
-
019
0.5000The instantaneous frequencies of all IMFs and OIMFs calculated using (5) were illustrated and compared in Figure 4. The results show that the instantaneous frequencies of all components are single-valued. In the instantaneous frequency distributions of IMFs, quite a few negative frequencies occur to the first IMF; very few negative frequencies exist in the fourth IMF; and no negative frequencies appear in the other IMF components. The negative frequencies lack of actual physical meanings, indicating that some deficiencies still remain when Hilbert transform is applied to the first IMF component, which will be further interpreted in Section 4. It could be concluded that, aside from the first IMF, other IMF components have good Hilbert property.Figure 4
Instantaneous frequencies of the IMFs and OIMFs.
(a)
(b)
(c)
(d)
(e)
(f)As we have known, thejth OIMF cj∗(t) is the linear combinations of the IMFs ci(t) (i=1,2,…,j). So whether an OIMF component bears good Hilbert property is influenced by the first several IMF components, causing that there appear lots of negative frequencies in the ninth OIMF. Another characteristic of EMD is the smoothness of the IMF components [29]. However, the orthogonalization makes the instantaneous frequencies of the high-order (low-frequency) OIMF components fluctuate quickly owing to the fact that they include a part of low-order (high-frequency) IMF components. So, it could be concluded that although the Gram-Schmidt method improves the orthogonality among the components, it brings some new difficulties. Both the smoothness and the Hilbert property of the components are influenced to some degree by such an orthogonalization process.In the above analyses, the orthogonalization was executed in accordance with the sequence from low-order components to high-order ones. So the low-order components degrade the good performance of high-order components, such as the smoothness, during the process of the orthogonalization. However, if the Gram-Schmidt orthogonalization is implemented with the opposite sequence, these problems will be resolved because the high-order IMF components with good smoothness and Hilbert property do not degrade the performance of other components. Figure5 presents the instantaneous frequency distributions of OIMFs obtained in this way. The results show that new OIMFs bear the same good performance with that of the IMFs. Of course, the deficiencies of the low-order components will still bring some difficulties of HHT in estimating the wind spectra, which will be discussed as an emphasis in Section 4.Figure 5
Instantaneous frequencies of the OIMFs.Figure5 shows that the instantaneous frequency of each OIMF of the wind data fluctuates around a certain mean value over time. Compared with the mean frequency, the transient frequencies of a low-order component (high-frequency component) are of relatively small amplitude. By contrast, the instantaneous frequencies of a low-frequency component vary over a relatively large amplitude range, indicating greater nonstationarity of such components.


## 4. Estimation of the Evolutionary Spectrum of a Typhoon Based on the HHT

### 4.1. Priestley Evolutionary Spectrum

Letx(t) denote a real nonstationary stochastic signal, which has a mean value of zero and a finite variance. At any arbitrary time, x(t) satisfies the following conditions: (20)Ext=0,Ext2<∞.The nonstationary signalx(t) can be expressed in spectral decomposition form as follows:(21)xt=∫-∞+∞At,ωeiωtdZω,where A(t,ω) is the time- and frequency-dependent modulating function. Z(ω) denotes a random process with orthogonal increments such that(22)EdZω=0,EdZω2=dμω.The autocorrelation function of the random processx(t) is(23)Rxxt,t+τ=Ex∗txt+τ=∫∞∞A∗t,ωAt+τ,ωeiωtdμω.Ifdμ(ω) is differentiable and dμ(ω)=S(ω)dω, then (23) can be written as(24)Rxxt,t+τ=Ex∗txt+τ=∫∞∞A∗t,ωAt+τ,ωeiωtSωdω.Whenτ=0, (24) becomes(25)Ex2t=∫-∞∞At,ω2Sωdω=∫-∞∞Sxxt,ωdω,where(26)Sxxt,ω=At,ω2Sω.Here, S(ω) is the power spectral density of the corresponding stationary process, Sxx(t,ω) is an evolutionary spectral density function or a time-varying power spectral density, and A(t,ω) is the modulation function. When A(t,ω) is constant, (26) will reduce to the spectral expression for a stationary signal.


### 4.2. Discrete Expression for the Estimation of the Evolutionary Spectrum Based on the HHT

The evolutionary spectrum of a nonstationary process is obtained by extending the concept of the spectral decomposition of a stationary random process. However, when the evolutionary spectrum is applied to describe an actual physical phenomenon, estimation of its parameters is not easy to achieve. Both the power spectral density and the frequency- and time-dependent modulating functions must be estimated. Based on the aforementioned OIMFs, one can deduce the following:(27)x2t=∑j=1ncj∗t2=∑j=1ncj∗t,ωt2=∑j=1ncj∗2t,ωt.In (27), x(t) is first expressed as the sum of the OIMF components. Given that the components cj∗(t) (j=1,2,…,n) possess the Hilbert property, they can be expanded in the time and frequency domains by applying the Hilbert transformation, thereby allowing the transient frequency ωt to be determined at any arbitrary time. Thus, the cj∗(t) (j=1,2,…,n) are written as cj∗(t,ωt) (j=1,2,…,n). The orthogonality property of the components is used in the last step of (27).In an actual engineering problem, one often obtains only one realization of a random processx(t). In such a case, the mean of x(t) is taken to be equal to its sample value; thus, one can deduce the following: (28)∫-∞∞Sxxt,ωdω=Ex2t=x2t=∑j=1ncj∗2t,ω.Given that negative frequencies are meaningless and the evolutionary spectrum is an even function, (28) can be further rewritten as(29)∫-∞∞Sxxt,ωdω=2∫0∞Sxxt,ωdω.At any instant of timet and any frequency ωk, one has(30)Sxxt,ωkΔω=12∑j=1ncj∗2t,ωk.Equations (29) and (30) indicate that the evolutionary spectrum of a nonstationary signal x(t) can be expressed in the form of an HHT spectrum, thereby resolving the difficulties in the choice of the modulation function for the estimation of the evolutionary spectrum.For a stationary or nonstationary signal, one can also define the time-averaged power spectrum as follows:(31)S-xxωk=12T∑i=1N∑j=1ncj∗2ti,ωkΔt.Equation (31) presents the expression for the time-averaged power spectrum based on the HHT, where T denotes the length of the entire time interval, N is the number of time points, and Δt is the time increment that is equal to the reciprocal of the sampling frequency.


### 4.3. Mode Mixing in the EMD of a Strong Typhoon

Here, the time-averaged power spectrum of the fluctuating wind data given in Figure1 is estimated using the HHT, and it is further compared with the Fourier power spectrum. The results are presented in Figure 6 and show that the spectra obtained using the two methods are consistent with each other in the low-frequency range. The Fourier spectrum exhibits a linear distribution in logarithmic coordinates, whereas the HHT spectrum decreases rapidly after 0.15 Hz. Regardless of the method used, the time-averaged energy or power spectrum should satisfy the law of conservation of energy. Because the spectra obtained using the two methods are nearly identical in the low-frequency range, they should also remain consistent in the high-frequency range. Therefore, in this case, it can be concluded that the HHT spectrum analysis fails to provide the true result in the high-frequency range. The method of HHT spectrum analysis still faces certain difficulties that must be overcome.Figure 6
Time-averaged power spectral density of the wind fluctuation data.As we know, a convergence criterion is required for a sifting process to stop. In the analysis presented above, thresholds ofθ1 = 0.05 and θ2 = 0.5 were used. It should be clarified whether a more rigorous criterion could improve the HHT spectrum at high frequencies. With the thresholds set to θ1 = 0.01 and θ2 = 0.1, the wind fluctuations data can be decomposed into 12 IMFs and one residue. However, the time-averaged power spectrum obtained from these IMFs remains nearly unchanged. Thus, a more rigorous criterion is found to play small role in improving the result of the HHT spectrum analysis.The upper frequency limit in Fourier spectrum analysis is equal to half the sampling frequencyfs. By analogy, we consider whether the maximum frequency of the HHT spectrum is also dependent on the sampling frequency of the signal and whether its upper limit is less than fs/2. In this paragraph, to investigate the influence of the sampling frequency on the upper frequency limit of the HHT spectrum, the time histories of turbulent wind data simulated using the harmonic synthesis method [30] are analysed. Suppose that the mean wind velocity is 30 m/s, the wind spectrum is the Kaimal spectrum [3], and the sampling frequencies used are 8 and 16 Hz. The simulated turbulent velocities are shown in Figure 7, and the HHT spectra of the simulated signals are provided in Figure 8. From Figure 8, it is evident that, at the same sampling frequency, the upper frequency limit of the HHT spectrum is considerably less than that of the Fourier spectrum. Although the use of a higher sampling frequency can increase the upper frequency limit of the HHT spectrum, this still fails to resolve the deficiency observed in the HHT spectral analysis.Figure 7
Simulated turbulent wind velocities for different sampling frequencies.
(a)
Sampling frequency = 8 Hz
(b)
Sampling frequency = 16 HzFigure 8
Power spectra of the simulated turbulent wind velocities.
(a)
Sampling frequency = 8 Hz
(b)
Sampling frequency = 16 HzAmong the IMFs or OIMFs of the wind data, the high-frequency components are of small amplitude. Thus, we consider whether some portion of these high-frequency, low-amplitude waves may be lost during the sifting processes because of the limitations of the precision of data storage in a computer. Thus, the Fourier transform was applied to the first four OIMFs to determine whether they still retain the high-frequency fluctuations. The results are shown in Figure9. The plots indicate that the first component is not a strict monocomponent, as expected, but rather still spans a broad frequency range. The second component also contains a small amount of high-frequency content that the Hilbert transform misses, whereas the time-averaged power spectra of the 3rd and 4th components from the Hilbert transform agree well with the corresponding Fourier spectra. Thus, the high-frequency fluctuations that we seek are hidden in the first two OIMFs and have not yet been extracted from them. Because of its own limitations, the Hilbert transform fails to produce the real spectra for these two OIMFs, particularly for the first OIMF, which also explains why a small number of negative frequencies appear in the first component in Figure 5. Thus far, we have demonstrated why the HHT fails to produce the true time-averaged power spectrum of a powerful typhoon, namely, because a mode-mixing phenomenon exists in the high-frequency IMFs or OIMFs of the wind data, which do not satisfy the “monocomponent” hypothesis.Figure 9
Time-averaged power spectral densities of the first four OIMFs of the wind fluctuation data.
(a)
First OIMF
(b)
Second OIMF
(c)
Third OIMF
(d)
Fourth OIMFThe phenomenon of mode mixing is one of the primary weaknesses of EMD. This phenomenon is defined as the case in which either a single IMF consists of signals at widely disparate scales or signals of a similar scale reside in different IMF components. To overcome the problem of mode mixing, an ensemble EMD (EEMD) method has been developed [31]. In this method, EMD is performed over an ensemble of the signal plus Gaussian white noise. The addition of white Gaussian noise can alleviate or completely resolve the mode-mixing problem. Torres et al. [32] proposed a complete EEMD (CEEMD) method with adaptive noise that can eliminate the residual noise from reconstructed signals in EEMD, thereby avoiding the possibility that different realizations of signal plus noise might yield different numbers of modes.The CEEMD method was applied to the turbulent wind data, and the Hilbert transform was then applied to obtain the power spectra of the complete ensemble IMF (CEIMF) components. In the EEMMD and CEEMD methods, the choices of the amplitude of the added noise and the number of ensemble trials affect the decomposition results. Wu and Huang [31] suggested that the amplitude of the added white noise be set as 0.2 times the standard deviation of the signal and that the number of trials be a few hundred. For the data dominated by low-frequency signals, larger amplitude of noise should be adopted [32]. Although other researches also discussed the choice of these parameters [33, 34], no general conclusions were provided. The suitable values still depend on the specific applications. In this work, the amplitude of the noise is set to 0.2 and 0.5 times the standard deviation of the wind fluctuation data, respectively. The corresponding noise standard deviations are 0.58 and 1.45. The number of trials is 500. The time-averaged power spectrum of the first CEIMF obtained using this method is presented in Figure 10. However, the results indicate that CEEMD still fails to eliminate the mode-mixing effect in the first IMF component of the wind data, and as a result, this method still cannot produce the true time-averaged power spectrum of the typhoon. The authors deem that although the EEMD and CEEMD methods are capable of solving the mode-mixing problem in intermittent signals and they cannot eliminate the mode mixing that arises in the case of small frequency ratios, thus, coping with mode mixing in EMD remains a challenge. To avoid the failure of the HHT to evaluate the true time-averaged power spectrum and evolutionary spectrum of a typhoon, a more effective method should be proposed, as discussed in the following sections.Figure 10
Time-averaged power spectral densities of the first CEIMF of the wind fluctuation data.
(a)
Noise standard deviationε=0.58
(b)
Noise standard deviationε=1.45


### 4.4. Analysis of the Time-Averaged Power Spectrum Based on the HHT-Fourier Transform

The degrees of nonstationarity of different-order OIMFs of wind data are different. Here, the variation coefficients of the variances in different time segments are used to compare the degrees of nonstationarity of different-order OIMFs, and the run test method described by Bendat et al. [35] is applied to determine whether a component could be deemed as stationary. We divided the wind fluctuation data in Figure 1(b) into 30 time segments; that is, each time of period of 2 minutes was treated as a short-time segment, and 30 variances were thus obtained. The degree of smoothness of these variances among the time segments reflects the degree of stationarity of the wind velocity. If the wind data were stationary, the variance would be constant. Table 2 provides the variation coefficients of the variances, which indicate that the long-cycle components exhibit strongly nonstationary features, whereas the high-frequency components exhibit good stability. Table 3 provides the round number [35] for each component and the acceptable range at a 10% level of significance. These results show that the assumption of stationarity is not rejected for the high-frequency components in the OIMFs but that some of the long-cycle components fail to pass the run test of stationarity.Table 2
Variation coefficients of the variances of all OIMFs.
Items
Variation coefficients
c
1
0.2601
c
2
0.5220
c
3
0.6193
c
4
0.9482
c
5
1.0616
c
6
1.0979
c
7
1.1892
c
8
1.3300
c
9
1.1070
r
9
0.8872Table 3
Run tests of the OIMFs.
Items
Round number
Acceptable range
c
1
15
[
11,20
]
c
2
18
[
11,20
]
c
3
16
[
11,20
]
c
4
14
[
11,20
]
c
5
15
[
11,20
]
c
6
9
[
11,20
]
c
7
11
[
11,20
]
c
8
8
[
11,20
]
c
9
11
[
11,20
]
r
9
15
[
11,20
]Based on the analysis of the nonstationarity degrees of the components, a comprehensive method of estimating the spectral distribution of a powerful typhoon is proposed. In this method, the low-frequency components with strong non-stationarity are analysed using the Hilbert transform, whereas the Fourier transform is applied to the high-frequency components that pass the run test of stationarity at a high level of significance. In this manner, we can retain the nonstationarity of a powerful typhoon but also avoid missing the high frequencies in the HHT time-averaged spectrum.For the data considered here, the Fourier transform was applied to the first two components and all other components were analysed using the Hilbert transform. After the time-averaged power spectra of the different-order components were estimated, the power spectrum of the original fluctuating wind can be obtained by superposing these spectra because of the complete orthogonality of all of the OIMF components. The analysis results, which are presented in Figure11, indicate that the trend of the time-averaged power spectrum is nearly identical to that of the Fourier spectrum. Thus, the time-averaged power spectrum of the typhoon has been successfully obtained, and the comparison with the Fourier spectrum indicates that it agrees well with the law of conservation of energy. The proposed method therefore resolves the deficiency of the HHT when used for spectral analysis due to the effect of mode mixing. The evolutionary spectrum of Typhoon Haikui will be presented in further detail below.Figure 11
Time-averaged power spectral density of the measured wind data.


### 4.5. Evolutionary Spectrum Analysis of a Strong Typhoon

In the previous section, the proposed method of spectral analysis based on the HHT was illustrated using one hour of wind data as an example. In this section, we study the evolutionary spectrum of the entire process of the typhoon, using all of the wind data from Typhoon Haikui as the object of analysis. Figure12 presents the overall wind data from Typhoon Haikui collected by the ultrasonic anemometer installed at the top of the south tower of the Sutong Yangtze Bridge. The mean velocity of the typhoon exhibits an overall trend of rising and falling. Here, however, we are interested in the time-varying power spectrum of the fluctuating wind during the whole process of the typhoon. Therefore, the mean velocity should be subtracted from the original wind velocity. The traditional method of obtaining wind fluctuations is to subtract the 10 min mean velocities from the wind data. However, this procedure will create some artificial “jumps” in the fluctuating wind data because different mean values are subtracted in different time segments [36]. Xu and Chen [37] proposed a method of extracting a continuously varying mean velocity to avoid such artificial “jumps” in fluctuating wind data. This method uses the residue of the IMF components to obtain the slowly varying mean. In this study, the method proposed by Xu and Chen [37] was applied to the wind data from Typhoon Haikui and the wind fluctuations were computed by subtracting the time-varying mean wind. The results of this calculation are presented in Figure 13.Figure 12
Wind velocity of “Haikui” at the Sutong Yangtze Bridge site.Figure 13
Mean wind and fluctuating wind of “Haikui” at the Sutong Yangtze Bridge site.
(a)
Mean wind
(b)
Fluctuating windFigure13 shows that the trend of the time-varying mean wind is well consistent with that of the 10 min mean velocities. However, the artificial “jumps” in the calculated wind fluctuations are eliminated, yielding a better representation of the turbulence characteristics of the typhoon. The figures show that even once the mean trend in the wind data from the typhoon has been eliminated, many abrupt instantaneous pulses are still observed in the strong typhoon, indicating that Typhoon Haikui exhibits strongly nonstationary characteristics. Thus, the evolutionary power spectrum of the typhoon should be investigated. The fluctuating wind data were decomposed into 17 IMF components and a residual component through EMD, and orthogonal components were then obtained using the Gram-Schmidt method. Next, the evolutionary spectrum of the wind data sequence was estimated using (30). The discrete Fourier transform was applied to the first three high-frequency components, whereas the other components were analysed using the HHT. Figure 14 presents the three-dimensional distribution of the time-varying power spectrum over the time-frequency domain. It can be observed that both the amplitude and frequency of this powerful typhoon vary over time, and thus, the power spectrum is time-varying. The power spectrum becomes stronger as the wind pulsation increases.Figure 14
Evolutionary power spectrum of Typhoon Haikui.Figure15 presents the distribution of the evolutionary spectrum in the time-frequency domain in the form of a contour map. The evolutionary spectrum presented in Figure 15 clearly illustrates the regularity of variation in the energy of the typhoon in both the time and frequency domains simultaneously. The wind power varies at each time and frequency point, reflecting the obvious nonstationarity of Typhoon Haikui. The traditional Fourier analysis cannot produce the time-varying power spectrum of such a typhoon. The figure clearly shows with each instantaneous increase in pulsation the power spectrum value also markedly increases and the energy distribution extends towards relatively high frequencies. During the first several hours represented in Figure 15, when the typhoon has not yet arrived at the bridge site, the fluctuation energy of the wind remains low; the turbulence energy then increases sharply with the arrival of the typhoon at the bridge site and remains at a high level until the typhoon has passed the bridge site, after which the turbulence energy recovers to a relatively low level. It can be concluded that the method proposed in this paper allows the nonstationarity of the turbulence in a typhoon to be effectively captured.Figure 15
Contour map of the evolutionary power spectrum of Typhoon Haikui (units: dB).


## 4.1. Priestley Evolutionary Spectrum

Letx(t) denote a real nonstationary stochastic signal, which has a mean value of zero and a finite variance. At any arbitrary time, x(t) satisfies the following conditions: (20)Ext=0,Ext2<∞.The nonstationary signalx(t) can be expressed in spectral decomposition form as follows:(21)xt=∫-∞+∞At,ωeiωtdZω,where A(t,ω) is the time- and frequency-dependent modulating function. Z(ω) denotes a random process with orthogonal increments such that(22)EdZω=0,EdZω2=dμω.The autocorrelation function of the random processx(t) is(23)Rxxt,t+τ=Ex∗txt+τ=∫∞∞A∗t,ωAt+τ,ωeiωtdμω.Ifdμ(ω) is differentiable and dμ(ω)=S(ω)dω, then (23) can be written as(24)Rxxt,t+τ=Ex∗txt+τ=∫∞∞A∗t,ωAt+τ,ωeiωtSωdω.Whenτ=0, (24) becomes(25)Ex2t=∫-∞∞At,ω2Sωdω=∫-∞∞Sxxt,ωdω,where(26)Sxxt,ω=At,ω2Sω.Here, S(ω) is the power spectral density of the corresponding stationary process, Sxx(t,ω) is an evolutionary spectral density function or a time-varying power spectral density, and A(t,ω) is the modulation function. When A(t,ω) is constant, (26) will reduce to the spectral expression for a stationary signal.


## 4.2. Discrete Expression for the Estimation of the Evolutionary Spectrum Based on the HHT

The evolutionary spectrum of a nonstationary process is obtained by extending the concept of the spectral decomposition of a stationary random process. However, when the evolutionary spectrum is applied to describe an actual physical phenomenon, estimation of its parameters is not easy to achieve. Both the power spectral density and the frequency- and time-dependent modulating functions must be estimated. Based on the aforementioned OIMFs, one can deduce the following:(27)x2t=∑j=1ncj∗t2=∑j=1ncj∗t,ωt2=∑j=1ncj∗2t,ωt.In (27), x(t) is first expressed as the sum of the OIMF components. Given that the components cj∗(t) (j=1,2,…,n) possess the Hilbert property, they can be expanded in the time and frequency domains by applying the Hilbert transformation, thereby allowing the transient frequency ωt to be determined at any arbitrary time. Thus, the cj∗(t) (j=1,2,…,n) are written as cj∗(t,ωt) (j=1,2,…,n). The orthogonality property of the components is used in the last step of (27).In an actual engineering problem, one often obtains only one realization of a random processx(t). In such a case, the mean of x(t) is taken to be equal to its sample value; thus, one can deduce the following: (28)∫-∞∞Sxxt,ωdω=Ex2t=x2t=∑j=1ncj∗2t,ω.Given that negative frequencies are meaningless and the evolutionary spectrum is an even function, (28) can be further rewritten as(29)∫-∞∞Sxxt,ωdω=2∫0∞Sxxt,ωdω.At any instant of timet and any frequency ωk, one has(30)Sxxt,ωkΔω=12∑j=1ncj∗2t,ωk.Equations (29) and (30) indicate that the evolutionary spectrum of a nonstationary signal x(t) can be expressed in the form of an HHT spectrum, thereby resolving the difficulties in the choice of the modulation function for the estimation of the evolutionary spectrum.For a stationary or nonstationary signal, one can also define the time-averaged power spectrum as follows:(31)S-xxωk=12T∑i=1N∑j=1ncj∗2ti,ωkΔt.Equation (31) presents the expression for the time-averaged power spectrum based on the HHT, where T denotes the length of the entire time interval, N is the number of time points, and Δt is the time increment that is equal to the reciprocal of the sampling frequency.


## 4.3. Mode Mixing in the EMD of a Strong Typhoon

Here, the time-averaged power spectrum of the fluctuating wind data given in Figure1 is estimated using the HHT, and it is further compared with the Fourier power spectrum. The results are presented in Figure 6 and show that the spectra obtained using the two methods are consistent with each other in the low-frequency range. The Fourier spectrum exhibits a linear distribution in logarithmic coordinates, whereas the HHT spectrum decreases rapidly after 0.15 Hz. Regardless of the method used, the time-averaged energy or power spectrum should satisfy the law of conservation of energy. Because the spectra obtained using the two methods are nearly identical in the low-frequency range, they should also remain consistent in the high-frequency range. Therefore, in this case, it can be concluded that the HHT spectrum analysis fails to provide the true result in the high-frequency range. The method of HHT spectrum analysis still faces certain difficulties that must be overcome.Figure 6
Time-averaged power spectral density of the wind fluctuation data.As we know, a convergence criterion is required for a sifting process to stop. In the analysis presented above, thresholds ofθ1 = 0.05 and θ2 = 0.5 were used. It should be clarified whether a more rigorous criterion could improve the HHT spectrum at high frequencies. With the thresholds set to θ1 = 0.01 and θ2 = 0.1, the wind fluctuations data can be decomposed into 12 IMFs and one residue. However, the time-averaged power spectrum obtained from these IMFs remains nearly unchanged. Thus, a more rigorous criterion is found to play small role in improving the result of the HHT spectrum analysis.The upper frequency limit in Fourier spectrum analysis is equal to half the sampling frequencyfs. By analogy, we consider whether the maximum frequency of the HHT spectrum is also dependent on the sampling frequency of the signal and whether its upper limit is less than fs/2. In this paragraph, to investigate the influence of the sampling frequency on the upper frequency limit of the HHT spectrum, the time histories of turbulent wind data simulated using the harmonic synthesis method [30] are analysed. Suppose that the mean wind velocity is 30 m/s, the wind spectrum is the Kaimal spectrum [3], and the sampling frequencies used are 8 and 16 Hz. The simulated turbulent velocities are shown in Figure 7, and the HHT spectra of the simulated signals are provided in Figure 8. From Figure 8, it is evident that, at the same sampling frequency, the upper frequency limit of the HHT spectrum is considerably less than that of the Fourier spectrum. Although the use of a higher sampling frequency can increase the upper frequency limit of the HHT spectrum, this still fails to resolve the deficiency observed in the HHT spectral analysis.Figure 7
Simulated turbulent wind velocities for different sampling frequencies.
(a)
Sampling frequency = 8 Hz
(b)
Sampling frequency = 16 HzFigure 8
Power spectra of the simulated turbulent wind velocities.
(a)
Sampling frequency = 8 Hz
(b)
Sampling frequency = 16 HzAmong the IMFs or OIMFs of the wind data, the high-frequency components are of small amplitude. Thus, we consider whether some portion of these high-frequency, low-amplitude waves may be lost during the sifting processes because of the limitations of the precision of data storage in a computer. Thus, the Fourier transform was applied to the first four OIMFs to determine whether they still retain the high-frequency fluctuations. The results are shown in Figure9. The plots indicate that the first component is not a strict monocomponent, as expected, but rather still spans a broad frequency range. The second component also contains a small amount of high-frequency content that the Hilbert transform misses, whereas the time-averaged power spectra of the 3rd and 4th components from the Hilbert transform agree well with the corresponding Fourier spectra. Thus, the high-frequency fluctuations that we seek are hidden in the first two OIMFs and have not yet been extracted from them. Because of its own limitations, the Hilbert transform fails to produce the real spectra for these two OIMFs, particularly for the first OIMF, which also explains why a small number of negative frequencies appear in the first component in Figure 5. Thus far, we have demonstrated why the HHT fails to produce the true time-averaged power spectrum of a powerful typhoon, namely, because a mode-mixing phenomenon exists in the high-frequency IMFs or OIMFs of the wind data, which do not satisfy the “monocomponent” hypothesis.Figure 9
Time-averaged power spectral densities of the first four OIMFs of the wind fluctuation data.
(a)
First OIMF
(b)
Second OIMF
(c)
Third OIMF
(d)
Fourth OIMFThe phenomenon of mode mixing is one of the primary weaknesses of EMD. This phenomenon is defined as the case in which either a single IMF consists of signals at widely disparate scales or signals of a similar scale reside in different IMF components. To overcome the problem of mode mixing, an ensemble EMD (EEMD) method has been developed [31]. In this method, EMD is performed over an ensemble of the signal plus Gaussian white noise. The addition of white Gaussian noise can alleviate or completely resolve the mode-mixing problem. Torres et al. [32] proposed a complete EEMD (CEEMD) method with adaptive noise that can eliminate the residual noise from reconstructed signals in EEMD, thereby avoiding the possibility that different realizations of signal plus noise might yield different numbers of modes.The CEEMD method was applied to the turbulent wind data, and the Hilbert transform was then applied to obtain the power spectra of the complete ensemble IMF (CEIMF) components. In the EEMMD and CEEMD methods, the choices of the amplitude of the added noise and the number of ensemble trials affect the decomposition results. Wu and Huang [31] suggested that the amplitude of the added white noise be set as 0.2 times the standard deviation of the signal and that the number of trials be a few hundred. For the data dominated by low-frequency signals, larger amplitude of noise should be adopted [32]. Although other researches also discussed the choice of these parameters [33, 34], no general conclusions were provided. The suitable values still depend on the specific applications. In this work, the amplitude of the noise is set to 0.2 and 0.5 times the standard deviation of the wind fluctuation data, respectively. The corresponding noise standard deviations are 0.58 and 1.45. The number of trials is 500. The time-averaged power spectrum of the first CEIMF obtained using this method is presented in Figure 10. However, the results indicate that CEEMD still fails to eliminate the mode-mixing effect in the first IMF component of the wind data, and as a result, this method still cannot produce the true time-averaged power spectrum of the typhoon. The authors deem that although the EEMD and CEEMD methods are capable of solving the mode-mixing problem in intermittent signals and they cannot eliminate the mode mixing that arises in the case of small frequency ratios, thus, coping with mode mixing in EMD remains a challenge. To avoid the failure of the HHT to evaluate the true time-averaged power spectrum and evolutionary spectrum of a typhoon, a more effective method should be proposed, as discussed in the following sections.Figure 10
Time-averaged power spectral densities of the first CEIMF of the wind fluctuation data.
(a)
Noise standard deviationε=0.58
(b)
Noise standard deviationε=1.45


## 4.4. Analysis of the Time-Averaged Power Spectrum Based on the HHT-Fourier Transform

The degrees of nonstationarity of different-order OIMFs of wind data are different. Here, the variation coefficients of the variances in different time segments are used to compare the degrees of nonstationarity of different-order OIMFs, and the run test method described by Bendat et al. [35] is applied to determine whether a component could be deemed as stationary. We divided the wind fluctuation data in Figure 1(b) into 30 time segments; that is, each time of period of 2 minutes was treated as a short-time segment, and 30 variances were thus obtained. The degree of smoothness of these variances among the time segments reflects the degree of stationarity of the wind velocity. If the wind data were stationary, the variance would be constant. Table 2 provides the variation coefficients of the variances, which indicate that the long-cycle components exhibit strongly nonstationary features, whereas the high-frequency components exhibit good stability. Table 3 provides the round number [35] for each component and the acceptable range at a 10% level of significance. These results show that the assumption of stationarity is not rejected for the high-frequency components in the OIMFs but that some of the long-cycle components fail to pass the run test of stationarity.Table 2
Variation coefficients of the variances of all OIMFs.
Items
Variation coefficients
c
1
0.2601
c
2
0.5220
c
3
0.6193
c
4
0.9482
c
5
1.0616
c
6
1.0979
c
7
1.1892
c
8
1.3300
c
9
1.1070
r
9
0.8872Table 3
Run tests of the OIMFs.
Items
Round number
Acceptable range
c
1
15
[
11,20
]
c
2
18
[
11,20
]
c
3
16
[
11,20
]
c
4
14
[
11,20
]
c
5
15
[
11,20
]
c
6
9
[
11,20
]
c
7
11
[
11,20
]
c
8
8
[
11,20
]
c
9
11
[
11,20
]
r
9
15
[
11,20
]Based on the analysis of the nonstationarity degrees of the components, a comprehensive method of estimating the spectral distribution of a powerful typhoon is proposed. In this method, the low-frequency components with strong non-stationarity are analysed using the Hilbert transform, whereas the Fourier transform is applied to the high-frequency components that pass the run test of stationarity at a high level of significance. In this manner, we can retain the nonstationarity of a powerful typhoon but also avoid missing the high frequencies in the HHT time-averaged spectrum.For the data considered here, the Fourier transform was applied to the first two components and all other components were analysed using the Hilbert transform. After the time-averaged power spectra of the different-order components were estimated, the power spectrum of the original fluctuating wind can be obtained by superposing these spectra because of the complete orthogonality of all of the OIMF components. The analysis results, which are presented in Figure11, indicate that the trend of the time-averaged power spectrum is nearly identical to that of the Fourier spectrum. Thus, the time-averaged power spectrum of the typhoon has been successfully obtained, and the comparison with the Fourier spectrum indicates that it agrees well with the law of conservation of energy. The proposed method therefore resolves the deficiency of the HHT when used for spectral analysis due to the effect of mode mixing. The evolutionary spectrum of Typhoon Haikui will be presented in further detail below.Figure 11
Time-averaged power spectral density of the measured wind data.


## 4.5. Evolutionary Spectrum Analysis of a Strong Typhoon

In the previous section, the proposed method of spectral analysis based on the HHT was illustrated using one hour of wind data as an example. In this section, we study the evolutionary spectrum of the entire process of the typhoon, using all of the wind data from Typhoon Haikui as the object of analysis. Figure12 presents the overall wind data from Typhoon Haikui collected by the ultrasonic anemometer installed at the top of the south tower of the Sutong Yangtze Bridge. The mean velocity of the typhoon exhibits an overall trend of rising and falling. Here, however, we are interested in the time-varying power spectrum of the fluctuating wind during the whole process of the typhoon. Therefore, the mean velocity should be subtracted from the original wind velocity. The traditional method of obtaining wind fluctuations is to subtract the 10 min mean velocities from the wind data. However, this procedure will create some artificial “jumps” in the fluctuating wind data because different mean values are subtracted in different time segments [36]. Xu and Chen [37] proposed a method of extracting a continuously varying mean velocity to avoid such artificial “jumps” in fluctuating wind data. This method uses the residue of the IMF components to obtain the slowly varying mean. In this study, the method proposed by Xu and Chen [37] was applied to the wind data from Typhoon Haikui and the wind fluctuations were computed by subtracting the time-varying mean wind. The results of this calculation are presented in Figure 13.Figure 12
Wind velocity of “Haikui” at the Sutong Yangtze Bridge site.Figure 13
Mean wind and fluctuating wind of “Haikui” at the Sutong Yangtze Bridge site.
(a)
Mean wind
(b)
Fluctuating windFigure13 shows that the trend of the time-varying mean wind is well consistent with that of the 10 min mean velocities. However, the artificial “jumps” in the calculated wind fluctuations are eliminated, yielding a better representation of the turbulence characteristics of the typhoon. The figures show that even once the mean trend in the wind data from the typhoon has been eliminated, many abrupt instantaneous pulses are still observed in the strong typhoon, indicating that Typhoon Haikui exhibits strongly nonstationary characteristics. Thus, the evolutionary power spectrum of the typhoon should be investigated. The fluctuating wind data were decomposed into 17 IMF components and a residual component through EMD, and orthogonal components were then obtained using the Gram-Schmidt method. Next, the evolutionary spectrum of the wind data sequence was estimated using (30). The discrete Fourier transform was applied to the first three high-frequency components, whereas the other components were analysed using the HHT. Figure 14 presents the three-dimensional distribution of the time-varying power spectrum over the time-frequency domain. It can be observed that both the amplitude and frequency of this powerful typhoon vary over time, and thus, the power spectrum is time-varying. The power spectrum becomes stronger as the wind pulsation increases.Figure 14
Evolutionary power spectrum of Typhoon Haikui.Figure15 presents the distribution of the evolutionary spectrum in the time-frequency domain in the form of a contour map. The evolutionary spectrum presented in Figure 15 clearly illustrates the regularity of variation in the energy of the typhoon in both the time and frequency domains simultaneously. The wind power varies at each time and frequency point, reflecting the obvious nonstationarity of Typhoon Haikui. The traditional Fourier analysis cannot produce the time-varying power spectrum of such a typhoon. The figure clearly shows with each instantaneous increase in pulsation the power spectrum value also markedly increases and the energy distribution extends towards relatively high frequencies. During the first several hours represented in Figure 15, when the typhoon has not yet arrived at the bridge site, the fluctuation energy of the wind remains low; the turbulence energy then increases sharply with the arrival of the typhoon at the bridge site and remains at a high level until the typhoon has passed the bridge site, after which the turbulence energy recovers to a relatively low level. It can be concluded that the method proposed in this paper allows the nonstationarity of the turbulence in a typhoon to be effectively captured.Figure 15
Contour map of the evolutionary power spectrum of Typhoon Haikui (units: dB).


## 5. Conclusions

This paper considered the evolutionary spectrum of Typhoon Haikui based on the HHT. The method of estimating the evolutionary spectrum using the HHT was rigorously derived. The mode-mixing effect that arises in the EMD of strong typhoon data was discussed in depth. Based on an analysis of the degrees of nonstationarity of the OIMFs, the paper proposed a synthesized method of studying the time-averaged spectrum and evolutionary spectrum of a strong typhoon:(1)
The IMF components that are obtained directly through EMD exhibit only approximate orthogonality. Energy leakage will occur if these IMF components are used to estimate the evolutionary spectrum and time-averaged power spectrum of a typhoon. Orthogonal components can be obtained by applying the Gram-Schmidt method to these IMFs. The sum of the squares of the OIMFs is consistent with that of the fluctuating wind velocity data, and no energy leakage occurs. During the process of the Gram-Schmidt orthogonalization, the first several components exert important influence on the performance, such as the smoothness and Hilbert property, of the orthogonalized components. Given that the high-order IMF components bear good performance, the orthogonalization should be implemented in accordance with the sequence from high-order components to low-order ones. In this way, the OIMFs bear the same good performance with that of the IMFs.(2)
The presented research proves that the EMD technique cannot effectively extract all high-frequency fluctuations from the wind velocity data of a strong typhoon; however, these waves remain in certain low-order (high-frequency) IMF components. Thus, these low-order components are not strictly “narrow-band” and do not possess a good Hilbert property. Therefore, the time-averaged power spectrum estimated using the HHT misses a portion of the high-frequency content of the data. To address this problem, this paper analysed the degrees of nonstationarity of different-order OIMFs and proposed a synthesized method of estimating the evolutionary spectrum and time-averaged spectrum. In this method, the Hilbert transform is applied to high-order components of strong nonstationarity and the Fourier transform is applied to low-order components that pass the run test of stationarity at a high level of significance. The time-averaged power spectrum obtained using this method is consistent with the Fourier spectrum.(3)
The evolutionary spectrum obtained based on the HHT can represent the energy of a powerful typhoon throughout the time-frequency domain and can clearly reflect the characteristics of the energy distribution over the time and frequency domains simultaneously. A large amount of instantaneously increasing pulsation is evident in the fluctuating velocity data from Typhoon Haikui, indicating significant nonstationarity. An analysis of the evolutionary spectrum indicates that when the pulsation of the typhoon increases, the evolutionary spectrum becomes stronger and the energy distribution spreads towards the relatively high-frequency range.This paper focused on providing an effective method of evaluating the evolutionary spectrum of a typhoon. However, only the longitudinal turbulent velocity was measured and analysed as an example of the application of the proposed method. In future research, a complete study of the turbulence characteristics of a typhoon, including both the longitudinal and vertical components, and the two-point coherence function of the different turbulence components will be performed. The method of eliminating the mode-mixing effect and abnormal instantaneous frequencies in the HHT is worthy of further study.



---

*Source: 4143715-2016-02-01.xml*