用于平稳过程的自适应多窗谱估计
工程技术TOPEI检索SCI升级版 工程技术1区IF 7.9高亮
- •
自适应多窗方法在频谱估计中采用每个频率点变化的窗数量。 - •
所提出的方法在两个具有不同光谱特性的例子中被证明是有效的。 - •
所提出的方法可以平衡光谱的局部分辨率和方差。
摘要
多窗谱估计方法在减少平稳随机过程谱的偏差和方差方面是有效的。该方法中的窗数量(NoT)控制着谱的偏差和方差之间的权衡。通常,NoT 是经验确定的,并且在所有频率下保持恒定,限制了局部偏差和方差的调整。本文提出了一种自适应多窗方法,该方法在每个频率点使用不同的 NoT 来估计多元平稳过程的功率谱密度(PSD)和相干函数。提出了一种带有停止准则的迭代过程,以在每个频率下优化 NoT,而无需手动调整。采用具有简单解析表达式和满意的泄漏保护的正弦窗。所提出的自适应多窗方法应用于三个数值示例:建筑物模型的结构响应、风速过程和自回归过程,这些示例表现出不同的谱特征。 所有示例的结果表明,所提出的方法在估计功率谱密度(PSD)和相干性时优于 Welch 方法和传统多谱图方法,具有更小的偏差和方差。
关键词
多谱图方法
功率谱密度
相干性
平稳过程
1. 引言
各种具有随机特征的物理和工程事件,如风速和在随机激励下的结构响应[1],[2],[3],[4],[5],可以视为平稳随机过程的一种实现。与时间域方法相比,通过频谱分析获得的频域表示更具信息量,允许分析能量分布[6],[7],[8]。
由于频域表示通常是从有限数据估计得出的,所获得的频谱通常会受到偏差和误差方差的影响[9]。几种非参数频谱估计器,包括 Welch 方法[10]和多谱估计方法[11],已被广泛用于减少估计频谱的偏差和方差[12],[13],[14],[15]。不幸的是,这些方法涉及偏差-方差权衡的困境,即低偏差估计器往往具有高方差,而低方差估计器通常具有高偏差[16]。多谱估计方法采用多个正交序列(加窗函数)进行频域平均以估计频谱[17]。 它相对于 Welch 方法表现出更好的性能,因为每个加窗的周期图都是从所有数据中获得的[18],[19],[20]。在实践中,Welch 方法可以被视为多窗谱估计技术的一种特殊情况,其中窗函数是相同的但时间上是移位的[20]。尽管这些方法可以缓解偏差和方差的问题,但它们通常包含难以调整的特定参数,例如 Welch 方法中的合适分段数和多窗方法中的窗函数数。
由自谱和互谱导出的相干性是量化两个时间序列在频域中线性耦合强度的重要工具[21]。频谱中的小低方差可能会导致较大的相干性方差。在工程应用中,准确估计相干性(如风相干性)具有挑战性,因为 Welch 方法和多谱估计方法是频率无关的,从而无法适应不同的频率分辨率[22]。
在最近几十年里,人们提出了几种方法来解决多谱估计法的偏差-方差权衡问题[19],[23],[24],[25],[26],[27],[28]。Thomson 提出了离散球面椭圆序列(也称为 Slepian 序列)作为一组具有较小泄漏的正交谱估计器,用于估计平稳过程的频谱[11],[29]。在设定时频积之后,可以计算出实现最佳方差的谱估计器数量。然而,控制频率分辨率和方差减少的带宽通常是由用户定义的,而没有系统的选取方法。 为解决这一问题,Haley 和 Anitescu [26] 提出了一种通过最小化对数谱的均方误差(MSE)来确定最优带宽的优化方法。尽管 Thomson 的多谱估计方法应用广泛,但 Slepian 窗的计算需要进行数值特征值分解,对于长时间序列来说计算成本较高 [30]。这种窗还需要在不同的带宽下重新计算。Riedel 和 Sidorenko [23] 提出了多谱估计中的正弦窗,它具有简单的解析表达式,无需设置和求解特征值方程。特别地,正弦窗的带宽仅取决于所用窗的数量。这意味着谱的分辨率和方差减少之间的内在权衡可以通过仅改变正弦窗的数量来控制。具体来说,较大的窗数量可以减少方差,但会伴随较高的偏差。 Riedel 和 Sidorenko [23] 还给出了通过最小化均方误差(MSE)来确定每个频率下使用的谱线数的公式。利用这个公式,Barbour 和 Parker [30] 开发了一种迭代程序来优化每个频率下采用的谱线数。尽管他们的方法生成了令人满意的频谱,但 迭代次数 是由用户根据经验确定的。最近,一种多谱线 S 变换(MTST)方法被提出用于 平稳和非平稳过程的谱估计,该方法结合了正交时频 Hermite 窗口和 S 变换方法 [27], [28]。该方法通过改变 Hermite 窗口的形状来降低不同频率下的谱和相干性的偏差和方差。 然而,MTST 方法在各种频谱形状上的适用性是有限的,因为其性能主要取决于定义的窗口形状函数,并且 NoT 在所有频率上仍然是恒定的。
在这项研究中,我们提出了一种自适应多谱估计方法,用于多变量平稳过程的谱估计和相干性估计,无需人工干预。与传统多谱估计方法固定 NoT 不同,我们提出了一种带有停止准则的优化策略来选择频率依赖的 NoT。因此,它使谱估计器能够在局部分辨率和谱方差之间取得平衡。通过三个数值例子验证了自适应多谱方法在处理不同谱特征方面的有效性,即建筑物模型的结构响应、风速过程和自回归(AR)过程。本方法减少了 NoT 选择的主观性,并平衡了局部谱的偏差-方差权衡。它还为具有各种谱形状的工程应用中的谱估计和相干性估计提供了一种有前景的方法。
2. 静态过程的多窗谱估计
2.1. 多变量平稳过程
设 为一个零均值多变量平稳过程,具有 , 。存在一个定义在区间 上的正交过程 ,使得[16] (1) 其中 T 表示转置算子, 为采样间隔, 为采样频率, 是奈奎斯特频率。 是一个复值正交增量过程,具有以下性质[16] (2) (3) (4) 其中*表示复共轭操作, 是克罗内克δ函数(当 时其值为 1,否则为 0),并且 被假定为纯粹连续的。
and the coherence function of is given by(6)where is the ijth element of .
2.2. Multi-taper spectral estimation
Tapering is an effective way to balance the broad- and narrow-band bias of the spectral estimation. Let be a discrete parameter stationary process with L observations. Assuming is the m-order taper, the multi-taper spectral estimator has the following form [27](7)and(8)where M is the maximum order of tapers. satisfies the orthogonal condition(9)
The selection of tapers significantly affects the resultant spectral estimation. The sine taper with a simple analytic expression is selected in this study. It is expressed as [23](10)
3. Adaptive multi-taper approach
3.1. Formula of the adaptive multi-taper method
As mentioned previously, the NoT of the multi-taper-based approaches controls the bias-variance trade-off in the spectral estimation. Usually, the NoT is fixed to be a constant integer empirically for all frequency points. The multi-taper spectral estimator with the frequency-dependent NoT is calculated as(11)
The estimator of the coherence is obtained as(12)
3.2. Bias and variance of the adaptive multi-taper method
with(14)where the spectral window has the requirement that .
The univariate stationary process is employed here for the convenience of demonstration. The bias of the multi-taper method is calculated by [16](15)
Assuming is double continuously differentiable and can be expanded in a second-order Taylor series with respect to f:(17)
of the sine taper is approximately rectangular over the resolution band , then [18], [26](18)where W is the bandwidth. Substituting Eqs. (17), (18) into Eq. (16), we can obtain(19)
If the spectrum is infinitely differentiable, the variance of the adaptive multi-taper method can be written as [18](20)
3.3. The optimal NoT
The MSE of the log spectrum estimated by the adaptive multi-taper method is given as [16], [26](21)where is used [26].
The optimal NoT at each frequency is determined by setting the differentiation of Eq. (21) with respect to M as zero:(22)
To obtain the stable , the ratio is transformed to the logarithm by setting(23)
Consequently, the first-order and second-order derivatives of y with respect to f are(24)(25)
Then, we can obtain the ratio as(26)
As a result, Eq. (22) becomes(27)
According to Eq. (27), to find the optimal , one needs to know the true log spectrum and its first- and second-order derivatives. An iterative procedure is proposed to obtain the optimal NoT as follows.
Step I: One taper for all frequencies is employed to generate the initial spectral estimate with a high resolution.
Step II: Fit a quadratic function to the estimated spectrum of the previous iteration within an interval of for each frequency f. The spectrum is extended by mirror symmetry to fit the spectrum near both ends of the frequency domain. The first- and second-order spectral derivatives of the log spectrum are then estimated and applied to Eq. (27) to calculate NoT.
Step III: From Eq. (22), it is noted that the NoT will increase to a large and even infinite value when the second-order derivative is small or vanish. In this regard, a forward and backward filtering technique is used to smooth the estimated NoT. A simple first-order low-pass filter is used as(28)where is the filter coefficient, is the original NoT before filtering at frequency f, and is the filtered value at f. The smaller the , the smoother the filtered spectrum is. In this study, is set as 0.1 to efficiently mitigate the distortion of the estimated NoT induced by the small second-order derivative of the spectrum. After that, the following constraint is placed on the filtered NoT to ensure the difference of the NoT at consecutive frequencies within one:(29)
Step IV: A two-step updating strategy is used to obtain the satisfactory NoT. First, repeat Steps II and III until the mean relative difference (MRD) of the estimated NoT between consecutive iterations is less than 10 %, or the number of iterations reaches 20. The MRD is defined as follows:(30)where J is the number of frequencies, and denotes the NoT at frequency of the sth iteration.
Step V: As fewer tapers are required for peaks in the spectrum to ensure a sufficient frequency resolution, the second updating step is adopted to reduce the NoT. In practice, the second-order derivatives of the spectrum peak are negatively large. Thus, the NoT is reduced by half at each iteration if the second-order derivative is less than , where and are the mean and standard deviation of the second-order derivative of the log spectrum, respectively. This correction step will be terminated until the MRD between consecutive iterations decreases to less than 1 %, or 30 iterations reach.
The flowchart of the method is displayed in Fig. 1.
The proposed method can calculate the optimal NoT at each frequency, which allows flexibility in adjusting the local bias and variance of various spectra. Three typical examples with different spectral characteristics are adopted to demonstrate the effectiveness of the method. The performance of the present adaptive multi-taper approach will be compared with Welch’s and traditional multi-taper methods.
4. Case study 1: A 4-DOF model
A four-story shear building model is used to evaluate the performance of the adaptive multi-taper method. The four degrees-of-freedom (DOF) structure is shown in Fig. 2, where the mass (ms) of each story is 1.0 ton, and the stiffness (ks) of each story is 180 kN/m. The natural frequencies of the building are 0.74, 2.13, 3.27, and 4.01 Hz.
The damping ratio of each mode is chosen to be 2 %. The building model is assumed to be subjected to a white noise excitation at each DOF. The spectral density matrix for the excitation is , where denotes a diagonal matrix. The responses of the model with a 400-s duration at a sampling frequency of 10 Hz are calculated by the 4-order Runge-Kutta method.
4.1. The spectrum estimate
The PSD of the structural responses is estimated using Welch’s, the traditional multi-taper, and the present adaptive multi-taper methods for comparison. The window length of Welch’s method and NoT of the traditional multi-taper method are determined such that both can provide their best results. Specifically, the window length of Welch’s method is 40 s with 50 % overlapping. The first ten sine tapers are used in the traditional multi-taper method. The NoT of the adaptive muti-taper approach is determined by the proposed optimization procedure. The optimization is terminated after 12 iterations, and the results are shown in Fig. 3. The NoT varies between 5 and 150. Fewer tapers are required at the PSD peaks, indicating a higher frequency resolution there.
The estimated PSDs of the structural responses at the first two floors are compared in Fig. 4. The Welch’s and multi-taper methods provide sufficient accuracy around the PSD peaks at the natural frequencies but have a large variance over the whole frequency range. The variances can be reduced by increasing the number of segments in Welch’s method and tapers in the traditional multi-taper method, which will increase the bias and decrease the frequency resolution. The proposed method leads to a similar resolution as the traditional multi-taper method at natural frequencies because it adopts a small NoT. It uses a larger NoT at other frequencies to reduce the variance.
The PSD estimates using the adaptive multi-taper method have smaller DSEs than those using Welch’s and traditional multi-taper methods over the whole frequency range.
The decibel-scale normalized relative error (DNRE) is then used to compare the errors quantitatively. It is expressed as [31](32)where and are the estimated and theoretical PSDs at fj, respectively, and denotes the Frobenius norm. In engineering applications, spectral moments are important to determine the survival probability and assess the reliability of structural systems [32]. The spectral moment is defined as [33](33)where p is the order of the spectral moment. The relative percentage error (RPE) between the estimated and theoretical spectral moments is given by(34)where is the estimated p-order spectral moment.
The DNRE, , and are listed in Table 1. The adaptive multi-taper method exhibits the smallest DNRE and RPE values in comparison with the other two methods. Hence, using the adaptive NoT results in more accurate PSD estimates than the traditional multi-taper method using a fixed NoT (ten tapers in this case) at all frequencies. The result also verifies the effectiveness of the proposed optimization strategy.
Empty Cell | Methods | DNRE (%) | (%) | (%) |
---|---|---|---|---|
First floor | Welch | 4.95 | 3.02 | 13.79 |
Traditional multi-taper | 6.29 | 2.37 | 12.68 | |
Adaptive multi-taper | 3.81 | 0.98 | 12.09 | |
Second floor | Welch | 4.64 | 4.77 | 4.88 |
Traditional multi-taper | 6.06 | 3.13 | 3.97 | |
Adaptive multi-taper | 3.70 | 2.47 | 3.38 |
4.2. The coherence estimate
The accuracy of the coherence estimate depends on the estimated PSD. Even a small fluctuation in the PSD could lead to a significant variance in the coherence estimate. The estimated coherence of the first two floors’ responses is presented in Fig. 6. The adaptive multi-taper method results in a more accurate coherence with a smaller fluctuation over the whole frequency domain than the other two methods. This is because increasing the NoT can also reduce the coherence variance, according to the statistical properties of the coherence of the multi-taper method [27]. In summary, the adaptive multi-taper method exhibits superior performance in balancing the bias and variance of the spectrum and coherence of structural responses.
4.3. Multiple realizations
Averaging multiple independent realizations can reduce the random fluctuation in the estimated spectrum and coherence. To demonstrate the stability of the optimal NoT and the effectiveness of the optimization procedure, the proposed method applies to 100 random realizations independently. According to the proposed adaptive approach, the optimal NoT for each realization is obtained. All 100 sets of the NoT are averaged and then compared with the theoretical one in Fig. 7. The theoretical NoT is calculated based on the theoretical PSD and Eq. (27). The figure shows that the averaged optimal NoT is close to the theoretical counterpart, except near the natural frequencies where the theoretical one dramatically increases to a high value. This is because the first-order and second-order derivatives of the PSD at these frequencies are approximately zero, as shown in Fig. 8. Accordingly, the optimal NoT becomes very large, according to Eq. (27). To avoid this happening, the difference of NoT at consecutive frequencies can be limited to less than one, as defined in Eq. (29). The constrained theoretical NoT is presented in Fig. 9.
The constrained theoretical and averaged optimal NoT are then used to estimate the PSD and coherence of the 100 samples. The averaged PSD and coherence over 100 realizations are calculated and presented in Fig. 10, Fig. 11. The 95 % confidence intervals of the results estimated by the adaptive multi-taper with the optimal NoT are also shown in these figures. Slight differences are observed in the estimated PSD using the theoretical and optimal NoT. The biases can be found between the averaged spectra and theoretical ones at the natural frequencies. This is because the adaptive multi-taper method minimizes the variance on the premise of sufficient resolution, as shown in Eq. (21). The confidence interval of the PSD and coherence estimates is narrow since the variance decreases with the increase of the NoT, according to Eq. (20).
The decibel-scale averaged absolute bias (AAB) and averaged coefficient of variation (ACoV) of the PSD and coherence, and RPE of spectral moments over 100 realizations are calculated and listed in Table 2. The AAB and ACoV are defined as(35)(36)
Estimate | Adaptive multi-taper | ||
---|---|---|---|
Theoretical NoT | Optimal NoT | ||
PSD for first floor | AAB | 0.358 | 0.359 |
ACoV | −0.078 | −0.088 | |
(%) | 4.425 | 4.816 | |
(%) | 0.915 | 1.954 | |
PSD for second floor | AAB | 0.353 | 0.359 |
ACoV | 0.027 | −0.004 | |
(%) | 5.776 | 5.828 | |
(%) | 2.939 | 3.736 | |
Coherence | AAB | 0.075 | 0.106 |
ACoV | −0.088 | −0.095 |
with(37)where denotes the kth realization at frequency , and K is the number of realizations (here is 100). AAB denotes the discrepancy between the mean of the estimate and the true value, while ACoV quantifies the averaged coefficient of variation. As displayed in Table 2, the AAB, ACoV, and RPE of the results estimated by the adaptive multi-taper method using the constrained theoretical NoT are close to those using the averaged optimal NoT. Given the above results, it can be concluded that the adaptive multi-taper approach can reduce the bias and variance of the PSD and coherence estimates.
5. Case study 2: A wind speed process
In this case, a bivariate stationary wind speed process is simulated and used to validate the adaptive multi-taper method. Its PSD is expressed as [27](38)