Characterizing non-Gaussian vibration loading using the trispectrum

This paper addresses the use of higher-order spectra to study the non-Gaussian nature of random vibration loading. Since the power spectral density is only a full description for stationary Gaussian processes, specifying non-Gaussian random vibration loading requires a sophisticated statistical description. In recent research higher-order statistical moments such as skewness and kurtosis have been used to define non-Gaussian properties of vibration loading. However, useful information contained in the spectral representation of these moments is neglected. This paper introduces the trispectrum as a tool for analyzing vibration loading. It is the spectral representation of the fourth-order moment and thus extends the information content of the kurtosis. For demonstration several common methods for generating non-Gaussian loading are reviewed and used to derive loads that reproduce the power spectral density and kurtosis of a real in-service loading. These loads are analyzed using Fatigue Damage Spectra and trispectra to relate structural response behavior to their non-Gaussian nature. The results suggest that the trispectrum is a valuable tool for analyzing and classifying non-Gaussian random loading.


Introduction
Designing reliable structures under the aspects of optimal dimensioning, a verified design and the saving of costs and time requires a comprehensive but efficient description of in-service loading. In numerous technical fields structures are subjected to random vibration loading, such as road, railway and aerospace vehicle excitation. Vibration loading must be analyzed within a fatigue strength assessment. Its purpose is to secure a structure's safe operation for a requested lifetime, making it an essential element throughout the design stages. Fatigue assessment under random vibration loading proves the difficulty of having a comprehensive load description available. This is due a loading's • random character.
• dependence on a set of constantly changing operating parameters (e.g. velocity, curves etc.).
• dependence on current environmental, structure and pavement conditions. Often there is no proper numeric model available that accounts for all variables and that is therefore capable of reproducing real loads adequately. Thus a common approach is to derive a statistical description from a realization of the vibration process of interest. This description may then be used to model the loading for a fatigue strength assessment of a new structure to be developed or for vibration qualification testing. If the in-service loading is of random nature it must be expressed by statistical methods. The derived statistical description must • reproduce structural responses that lead to the relevant failure mechanisms of a structure.
• be representative for all relevant characteristics of in-service loading. • be efficient to be assessed in a stress response analysis within a reasonable period of time.
• assign for the randomness. A quantity effectively incorporating these requirements is the power spectral density (PSD). The PSD depicts the amplitudes of all harmonics contained in a load series and its integral provides an estimate of a loading's variance. For Gaussian stationary random processes the PSD covers a full stochastic description. Measured in-field vehicle vibration usually differ significantly from this assumption. They are subjected to a variety of superimposed effects which interfere with stationary Gaussian behavior. In particular the constant change of conditions and parameters conflicts with the stationary hypothesis. The PSD of such complex loading determines the average power spectral density. A replication of the loading solely based on the averaged PSD generally has lower amplitudes than there are contained in measured load series. This results in non-conservative deviations for the lifetime estimation. As a consequence different techniques have been developed that aim to adequately abstract the relevant characteristics of random vibration loading [1]. Most of them extend the description of the PSD by additional statistical quantities. The central idea is to reproduce these measures in a synthetic load in order to subsequently reproduce real failure mechanisms in a fatigue assessment. In recent research this has primarily been the kurtosis. The kurtosis is a normalized higher-order moment that follows the variance in describing symmetric spread of random processes. Nevertheless statistical load description applies to a structure's excitation. But fatigue load spectra are derived from a structure's response and little attention has been payed to the transfer behavior of the non-Gaussian load character. In [2] it was observed that a structure's response subjected to non-Gaussian loading depends decisively on the non-Gaussian nature. Simple structures were simulated and tested with different synthetic non-Gaussian loads having the same PSD and similar kurtosis values [3,4]. The authors assessed strongly differing fatigue behavior. Others find in random loading deviations for the time-to-failure depending on the frequency bands adjusted to generate a specific kurtosis [5]. These findings indicate that the scalar value of the kurtosis is insufficient to adequately describe non-Gaussian random vibration loading. This paper pursues the idea that it is desirable to extend the load description by the PSD with a spectral representation of the non-Gaussian character. Higher-order spectra (HOS) lay the foundation for decomposing higher-order moments in frequency domain. The use of the 4 thorder spectrum, the trispectrum, is based on the motivation to provide a spectral representation of the kurtosis. It extends the non-Gaussian information of a loading and allows to identify and determine underlying mechanisms in detail. This paper starts with an introduction to statistical load analysis. From well-known higher-order quantities it is moved onto HOS and the trispectrum, associated with the spectral decomposition of the kurtosis. The following section reviews a set of methods to generate non-Gaussian loading. They are applied to derive synthetic loads according to the PSD and kurtosis of a real load series. The Fatigue Damage Spectrum (FDS) is used to calculate the response behavior of all loads using one-degree-of-freedom dynamic systems. These results are compared and related to the trispectra.

Statistical load analysis
Theoretically any stochastic process X(t) can be completely described by the full set of probability density functions p(x, t) (PDF). However, this goes with an effort and a complexity that cannot be realized in practical applications. One comes closer to practice by the assumption of stationarity. This implies that the PDF p(x) does not change when shifted in time. Nevertheless stationarity statistical quantities are derived by ensemble averaging. This requires a large amount of realizations, an ensemble, which are averaged to abstract the quantities of interest. Assuming ergodicity allows to derive these quantities from a single realization, by time averaging. This assumption lays the basis for practical statistical load descriptions. Ergodic RASD IOP Conf. Series: Journal of Physics: Conf. Series 1264 (2019) 012040 IOP Publishing doi:10.1088/1742-6596/1264/1/012040 3 processes can be described by a single PDF p(x) independently of time. A process is considered to be a random Gaussian process if it follows the parametric PDF p g (x).
Gaussian processes are fully specified by mean µ x and variance σ 2 . These are representatives of a more general class of statistical quantities -the statistical moments. Moments m n include the mean µ x = m 1 , while central moments µ n are measures about the mean.
The variance σ 2 = µ 2 is the basic measure for spread and for that reason is also used to normalize higher-order (order n > 2) central moments. The normalized moments following are the 3 rd -order skewness γ and the 4 th -order kurtosis β.
Higher-order central moments of even order are measures of symmetric spread, while unevens describe asymmetric spread. In theory an infinite set of statistical moments also defines a PDF. In practice, statistical moments can be used to identify or to quantify deviations from a parametric PDF. If a zero-mean process follows the Gaussian distribution its central moments are given in terms of µ 2 .
Thus cumulants c n are introduced. Cumulants are an alternative to moments that in some theoretical problems offer a favorable mathematical approach. Here they are presented due to their properties for the Gaussian distribution, where in contrast to central moments µ n (Eq. 4) all Gaussian cumulants of higher-order are zero c n = 0; n ∈ {3, 4, ..}. They provide a clear picture of the deviation to a Gaussian distribution and also to cumulants of other order. The normalized 4 th -order cumulant c 4 is more commonly known as excess kurtosis. The power spectral density (PSD) S xx (f ) is a fundamental tool of signal processing. Above all it might be the most essential tool in vibration fatigue. This is due to its • efficient load description by a single continuous function that e.g. allows it to be used for test standards. • clear physical meaning that enables to forecast or to analyze the vibration response of mechanical structures. • efficient processing in a stress response analysis of linear structures that allows to realize fast simulations in all design stages. • scalability and opportunity to derive an infinite amount of Gaussian realizations to carry out statistically robust analyzes. • possibility to estimate statistical distributions of load spectra for Gaussian responses. The PSD is related to the 2 nd -order central moment via Parseval's theorem. The integral of the PSD equals the 2 nd -order central moment. Therefore the PSD can be interpreted as the spectral decomposition of µ 2 . It indicates how each frequency contributes to the variance. This leads to the central idea of the following section, which pursues the idea to assign a spectral decomposition for the higher-order moments that are needed to describe non-Gaussian processes.

Higher-order spectra (HOS)
From the central limit theorem in probabilty theory follows that the Gaussian distribution is the general result of a large number of random process variables acting independently together. By contrast many real processes are influenced by a set of process mechanisms that interact to a greater or lesser extent usually hidden to the observer. In time domain these interactions may be assessed by auto-correlation functions R (x) n (τ 1 , ..., τ n−1 ). An equivalent description in frequency-domain is the spectral analysis S (x) n (f 1 , ..., f n−1 ), depicting the associated frequencies [6]. The spectral analysis generally processes the information of the Fourier spectrum X(f ).
It is a well-known fact that the 2 nd -order spectrum, the PSD S .} can be understood as extensions of the PSD, but include the phase information. They depict a time series not only by contributions of each single frequency, but also by cross-correlations between frequencies. Thus HOS have multiple arguments, indicating the cross-correlation structure in frequency domain. The n th -order spectrum refers to n−1 frequency arguments. Integrating this frequency space ends up to the corresponding higher-order central moment. This can be understood as an extension of Parseval's theorem (Eq. 5). (7) HOS decompose higher-order moments in frequency domain. This allows to give popular measures such as skewness and kurtosis a spectral decomposition. HOS extend their scalar value to a spectrum of frequency cross-correlations. By definition stationary Gaussian processes have no such correlation structure. HOS of Gaussian processes contain no additional contributions. The first HOS is the bispectrum S . The appendix -bi denominates the decomposition into a two-dimensional space, indicating frequency cross-correlations that form the 3 rd -order central moment. Normalizing S xxx (f 1 , f 2 ) leads to the spectral representation of the skewness. Symmetric processes have a skewness of γ = 0. The trispectrum S represents the decomposition of the 4 th -order central moment over the three-dimensional frequency space. Integrating a normalized trispectrum results in the kurtosis.
is the decomposition of the fourth cumulant c 4 or excess kurtosis. Bi-and trispectrum have been the most popular HOS as being the consecutive spectra to the PSD and being associated with principal asymmetry and spreading properties. HOS may be used for a detailed analysis of the frequency-domain content to allow to classify and quantify random processes and to study nonlinear systems [7]. Due to its potential there has been extensive research on the use of HOS [8]. However, few practical applications issued from these activities. This certainly can be related to the fact that the complexity of HOS goes far beyond the PSD. Redundancies, complex-values and multidimensionality make the interpretation, visualization and estimation of HOS challenging.

Non-Gaussian vibration loading
In the prior section a variety of stationary quantities are introduced that are generally estimated employing the assumption of ergodicity. This is a dilemma because measured realizations of vehicle vibration loading usually do not conform with the strict mathematical requirements. The result is a conflict between the strict mathematical formulation and the efficient measures of the reviewed stationary quantities. At the core is that the PSD is only a complete descriptor for stationary Gaussian loading. Techniques that have been developed to deal with non-Gaussian loading approach from two different points of view, either dropping one of the PSD's postulations. They can therefore be distinguished into non-stationary and non-Gaussian approaches. The former is based on the assumption that a real load series is generated by one resp. a set of Gaussian processes that vary in time. The corresponding methods aim to decompose the non-stationary load series into different stationary Gaussian segments. This semi-stationary approach allows to tie in the very efficient frequency methods (e.g. Dirlik-method) for fatigue assessment [9,10]. Nevertheless decomposing the load is challenging because the variation of conditions occur within one realization of the process. In contrast the non-Gaussian approach is of simpler use. The underlying assumption is that varying load levels occur frequently. So a load series is considered to be a stationary realization of the random process in question if it is observed until all its varying characteristics repeat themselves. Ergo it is chosen to be representative for the in-service operation. The consequential statistical description is based on stationarity but a non-Gaussian distribution. To define the non-Gaussian nature the description is extended to higher-order statistical quantities. The following review of methods are based on this non-Gaussian approach, reproducing the kurtosis. In this context the trispectrum is used to abstract more detailed information about the non-Gaussian character. It reveals differences and similarities between the load series. Knowing how a load series is composed may allow to describe it very efficiently by stationary quantities. Despite the differing handling of stationarity and their conflict with its strict mathematical formulation both approaches may be legitimated for fatigue assessment in the high-cycle range. This is due to the subsequent use of load cycle counting schemes (e.g. rainflow-counting) which derive load spectra independently of stationarity and load sequence effects. Synthetic non-Gaussian loading All methods to be introduced that generate non-Gaussian vibration loading begin with a stationary realization of a Gaussian process. This is the easiest way to generate a loading with a specific frequency content. The Fourier amplitudes are derived from a PSD and paired with random phases, statistically independent and uniformly distributed. This serves as the basis for carrying out manipulations of different kinds to make the loads non-Gaussian. The following methods are based on manipulating the phases, functional transforming the Gaussian PDF or modulating of loads.

Nonlinear Transformation via Hermite polynomials (HP)
In [11] Winterstein proposes a transformation G{·} based on Hermite polynomials (HP) that generates a process Y (t) with a nominal kurtosis β y . To adjust the kurtosis it is sufficient to use the third-order polynomial and determine its pre-factorh 4 (Eq. 8). In its basic set-up the transformation requires a normalized Gaussian input X(t) ∼ N (µ = 0, σ 2 = 1). To further handle the output standard deviation, the correction factor κ is needed for equally producing a normalized output.h Then the output can be linearly scaled to the desired standard deviation σ y . This generates a transformed process Y (t) with a nominal kurtosis β y and standard deviation σ y from a general zero-mean Gaussian input X(t).
The HP transformation is of simple use and the kurtosis serves as the single parameter to be adjusted. Nevertheless in the following the trispectrum is used to show that the transformation produces super harmonics that may interfere with the uniqueness of the PSD (sampling theorem).

Implementation for this paper
The adjustment of the kurtosis by formulation (Eq. 9) is an approximation that drifts with increasing kurtosis. It requires some iterations to approach the desired kurtosis.

Analytical Phase Selection (APS)
As discussed before, the PSD contains the full information of the Fourier-amplitudes, but is independent of the phases. In contrast HOS depend on the phase information. The idea of Analytical Phase Selection (APS) is in this context very elegant. The PSD is kept untouched while a manipulation of the phases is used to increase the fourth-order moment and thus the kurtosis (compare Eq. 3). Stationary Gaussian load series have a statistically independent and uniformly distributed phase structure. Via APS one tries to change portions of this structure to be deterministic in order to produce a desired kurtosis. Early propositions [12] were based on very simple manipulations, such as turning random phases into zero-phases until a desired kurtosis is reached. This increase of kurtosis is based on the fact, that harmonics with same phases add up on specific occasions. Nevertheless such a procedure generates load series that are far from reality. In [13] Steinwolf presents a more sophisticated method to generate phase dependencies. This method introduces deterministic phases that complement random phases to generate a pseudo-random phase structure. The theoretical framework of the method is strong but the manipulations presented in the publication are either not effective ('phase quartets') or its implementation is complex ('b k -method'). Partly because the explanation are not thorough to be realized in application. Still the concept of using phase manipulations to increase higher-order moments independently of the PSD is very appealing. It is desirable that a general applicable framework would be developed. Implementation for this paper To implement the APS method, the b k -method is used according to [13] with 50 phase combinations for each deterministic selected phase. In order to effectively generate a relevant kurtosis it was necessary to divide the load series ( Fig. 1) into blocks. The block size was chosen to be twice the sampling rate f s , equaling a time of two seconds. Nevertheless this segmentation leads to discontinuities and thus distortions in the PSD. To reduce this effect it needs a subsequent overlapping scheme. The phase manipulation also had to be performed on previously manipulated phases. To reduce the method to few parameters, the manipulation was carried out along the full frequency axis. Via the step size of deterministic chosen phases it was iterated until a specific kurtosis was reached.

Carrier-Wave (CW)
The principal idea of the carrier-wave method (CW) is that in reality vibration processes often vary in intensity due to changing operating, environmental or excitation source conditions. The CW method functions via amplitude-modulation, where a generally a low-frequency signal modulates a high-frequency carrier wave. For the purpose of generating a varying load the

Discussion
For an analysis of the introduced methods, realizations of load series are generated according to the prior section that reproduce PSD and kurtosis of a real in-service load (Fig. 1). This reference load was measured as an acceleration in vertical direction on the axle box of a railway locomotive. In comparison, the real load is characterized by the most irregular progression and a small positive skewness γ, common for railway loading. For all non-Gaussian methods, the kurtosis β is successfully adjusted to a value close to the the reference load. Also the graphs of the PSD (Welch-method, ∆f = 1 Hz) show conformity. Only the load generated by APS stands out with a more fluctuating PSD, which can be attributed to the block-approach. To compare the loads in terms of structural responses the FDS is used (Fig. 2). It simulates a general structural response behavior by varying the natural frequency f D of a one-degree-of-freedom system [14]. In addition, the Miner exponent k is varied to represent different structural fatigue behavior. The damping ratio D is kept constant. Every data point on the FDS functions represents a unique response analysis. Each producing a load spectrum with cycle amplitudes y j and the corresponding counts N j . For the FDS these are summarized as an equivalent fatigue damage d eq using the Basquin equation for a single cycle (Eq. 10). The structural response is proportional to this fatigue damage when applying Miner's rule.
As expected the Gaussian load series produces generally the least damage. The HP load series causes slightly more damage and stays relatively constant over frequency. The FDS of the loads generated via APS and CW are very similar. The small deviations may result from the fluctuating PSD of the APS load. Both follow the shape of the Gaussian FDS-curves but are shifted upwards. They scale the damage upwards by a factor that is also influenced by the Miner exponent. The measured load series generally produces the most damage, which appears more pronounced for larger Miner exponents. Especially the interval until 50 Hz and the intervals around the peaks at 160 Hz and 320 Hz stick out. For structures having natural frequencies in these intervals, large deviations in fatigue damage can be prognosticated, when replacing the real by the synthetic loads. Then the synthetic loads lack to reproduce the relevant failures mechanisms in a comparable time frame. This is very unsatisfactory, since neither the large deviations could have been foreseen from the equal kurtosis values nor the addressed intervals could have been expected on the basis of the PSD. The illustration of the trispectra is more difficult, since they include tremendously more data points that distribute on three dimensions. Figure 3(a) shows a qualitative plot of the cumulant trispectrum of the measured load series. This three-dimensional surface plot shows all contributions to the excess kurtosis, weighted in size. The cumulant trispectrum is chosen as it features deviations from the Gaussian distribution. Thus the Gaussian load series does not need to be included in the plots, being zero for the full frequency space. Nevertheless the three-dimensional illustration is very difficult to interpret. Thus for an analysis slices are used, of which two are discussed in the following. The locations of the slices are illustrated in Figure 3(b) and are taken along the f 3 -axis. Fig. 4 and 5 display these slices for all non-Gaussian loads. Every pixel represents an averaged cube of the trispectral space with an edge length of 5 Hz. All symmetries of the trispectra are included. The slices are chosen to highlight specific information of the trispectra.  The shape of the PSD function can be detected along these axes in grey tones -indicating a strong correlation between the 2 nd -and 4 th -order spectral densities. The trispectral slice of the measured load series extends along the same axes but is more centered around the origin. Hypothetically summing up the discussed slices, the measured load series would have a lower contribution than the other two loads. Since all loads have the same kurtosis the difference must be distributed on different slices. Consequently the trispectrum of the measured load provides new insight. The 4 th -order spectral density is distributed non-proportional to the PSD. In contrast the HP trispectral slice shows a geometric pattern. The intersections of the three axes form a triangle shape that is surrounded by a reversed triangle shape having the triple distance from the origin. This can be traced back to the negative term in Eq. 9. In this context the trispectrum may serve as a tool for system analysis of the transformation. In terms of understanding the differing structural response behavior it is more interesting to study the second slice in Figure 5. The straightforward conclusion is that only the real load series shows a considerable contribution to the excess kurtosis c 4 that can be attributed to the interval ∆f 3 = 155−160 Hz. A fact that correlates with the FDS (Fig. 2). Also the contributions at f 1 , f 2 ≈ 320 Hz need to be mentioned. They correlate with the second peak in the FDS. Due to the symmetries, also the slices corresponding to the second peak have these considerable contributions, but are not plotted here. Most of all, this finding implies a strong correlation between the two intervals. This suggest that the FDS, which operates on the single natural frequency of a one-degree-of-freedom system, may not be suitable to fully assess the response behavior of real structures. Real structures can have multiple resonance frequencies, which may be sensitive to this coupling. Such a coupling is not covered by the FDS. This section showed that non-Gaussian structural behavior can be related to the trispectrum. Another important aspect is shown in the lower slices of Fig. 5. In these the upper right corner is highlighted with a lower scale. It is expected that these sections do not include any relevant contributions since this would require a coupling with a fourth frequency f 4 = |f 1 + f 2 + f 3 | (Eq. 6) that is larger than the Nyquist frequency f 4 > f N y = fs 2 = 600 Hz (compare [15]). Nevertheless in the HP slice diagonal contributions corresponding to frequencies f 4 ≈ 1150 Hz can be identified and thus suggest the existence of higher frequency content. This is the result of the nonlinear transformation (super harmonics) and cannot be uniquely assessed by the PSD. But it represents a plausible explanation for the small damage caused by the HP load -a significant contribution is distributed on higher frequencies which are filtered out in the response analysis. Looking back at the far right corner of the FDS the HP load just starts surpassing the other loads, corroborating this assumption.

Conclusion
This paper addresses the use of higher-order spectra to study non-Gaussian loading and its effects on structural responses. This is demonstrated by analyzing different non-Gaussian loads using Fatigue Damage Spectra and trispectra. All loads are based on the same PSD and kurtosis that originated from a measured in-service load. Despite their conforming kurtosis values the trispectrum reveals different decompositions of the kurtosis. While all synthetic loads have trispectra that are proportional to the PSD shape, the real load series shows a trispectrum with non-proportional deviations. The corresponding frequency intervals are retrieved in the FDS, clearly stating a correlation between the trispectrum and structural responses. These findings show that the trispectrum is entitled as an instrument for the fatigue assessment of non-Gaussian random loading. Its application allows a wide field of other usages, which is exemplified by the identification of aliasing in the HP method. Also its application leaves numerous open questions. Two of them are how the correlation of frequencies affect the response of more complex mechanical systems and how the trispectrum can be reproduced in synthetic loads.