振动工程学报

非平稳非高斯随机过程插值模拟方法

PDF下载

盛向前 ¹ , 虞跨海 ¹ , 范文亮 ² , 牛兰杰 ³

振动工程学报 | 2025,38(8): 1827-1838

收起

振动工程学报 | 2025, 38(8): 1827-1838

非平稳非高斯随机过程插值模拟方法

全屏

盛向前¹, 虞跨海¹, 范文亮², 牛兰杰³

作者信息

^1.河南科技大学工程力学系，河南洛阳 471000

^2.重庆大学建筑力学系，重庆 400045

^3.西安机电信息研究所，陕西西安 710065

盛向前（1988—），男，博士，讲师。E-mail: 2008shengxq@163.com

通讯作者:

虞跨海（1982—），男，博士，教授。E-mail: yukuahai@163.com

Interpolation simulation method of non-stationary non-Gaussian stochastic processes

Xiangqian SHENG¹, Kuahai YU¹, Wenliang FAN², Lanjie NIU³

Affiliations

^1.Department of Engineering Mechanics，Henan University of Science and Technology，Luoyang 471000，China

^2.Department of Architectural Mechanics，Chongqing University，Chongqing 400045，China

^3.Xi’an Institute of Electromechanical Information Technology，Xi’an 710065，China

出版时间: 2025-08-10 doi: 10.16385/j.cnki.issn.1004-4523.202309035

文章导航

摘要

收起

针对非平稳非高斯随机过程模拟中存在的随机变量数目过多和潜在高斯随机过程的功率谱计算耗时大的问题，本文结合随机谐和函数，提出一种基于样本插值的非平稳非高斯随机过程快速模拟方法。在已知非高斯随机过程的目标演变功率谱和目标密度函数的前提下，通过Mehler公式建立非高斯随机过程和潜在高斯随机过程的相关函数方程，并采用插值求解的方式提出潜在高斯随机过程的演变功率谱快速计算方法，结合随机谐和函数提出非平稳非高斯随机过程快速模拟方法，采用单点均匀调制非高斯随机过程和多点非均匀调制非高斯随机过程模拟验证该方法的有效性。结果表明：在保证计算精度的前提下，插值求解潜在高斯随机过程的演变功率谱的计算耗时低于Mehler公式求解的耗时，且随着激励数目的增多，插值求解计算潜在高斯随机过程的演变功率谱的效率更为明显；所提非平稳非高斯随机过程快速计算方法，能够有效模拟具有目标演变功率谱和目标密度函数的非高斯随机过程。

关键词

随机振动 / 非高斯随机过程 / 非平稳 / 演变功率谱 / 随机谐和函数

Abstract

收起

To address the problems of large number of random variables and time-consuming computation in the simulation of non-stationary non-Gaussian stochastic processes，a fast computation method of non-stationary non-Gaussian stochastic processes is proposed based on sample interpolation by combining the stochastic harmonic function. With the known of the target evolutionary power spectrum and target density function of non-Gaussian stochastic processes，the correlation function equations of non-Gaussian stochastic processes and underlying Gaussian stochastic processes are established through Mehler’s formula，and a fast calculation method for the evolutionary power spectrum of underlying Gaussian stochastic processes is proposed through interpolation method. Subsequently，a fast simulation method for non-stationary non-Gaussian stochastic processes is proposed by combining stochastic harmonic functions，The effectiveness of this method is verified by simulating single-point uniformly modulated non-Gaussian stochastic process and multi-point non-uniformly modulated non-Gaussian stochastic processes. The results show that，when calculating the evolutionary power spectrum of the underlying Gaussian random process under the condition of ensuring accuracy，the calculation time of interpolation solution is lower than that of Mehler’s formula solution，and as the number of excitations increases，the efficiency of interpolation solution in calculating the evolutionary power spectrum of the underlying Gaussian random process is more obvious. The proposed fast computational method of non-stationary non-Gaussian stochastic processes can effectively simulate the non-Gaussian stochastic processes with the target evolutionary power spectrum and the target density function.

Key words

random vibration / non-Gaussian stochastic processes / non-stationary / evolutionary power spectrum / stochastic harmonic functions

引用本文

盛向前, 虞跨海, 范文亮, 牛兰杰. 非平稳非高斯随机过程插值模拟方法. 振动工程学报, 2025 , 38 (8) : 1827 -1838 . DOI: 10.16385/j.cnki.issn.1004-4523.202309035

Xiangqian SHENG, Kuahai YU, Wenliang FAN, Lanjie NIU. Interpolation simulation method of non-stationary non-Gaussian stochastic processes[J]. Journal of Vibration Engineering, 2025 , 38 (8) : 1827 -1838 . DOI: 10.16385/j.cnki.issn.1004-4523.202309035

正文

收起

工程结构在服役期间所遭受的荷载具有非高斯特性的现象，如地震^[1]、风压^[2]和火箭起飞的振动环境^[3]等逐渐引起研究者关注。为评估工程结构在非高斯随机激励下的响应特性，通常需要采用时域方法获得结构响应的时程信息。因此，在随机振动分析中，非高斯随机过程模拟尤为重要。

针对非高斯随机过程模拟，目前已发展多种方法，主要包括传递函数的谱表示法^[4-12]、Karhunen-Loeve （K-L）展开法^[13]和高阶谱表示法^[14]等，其中传递函数的谱表示法应用最为广泛。关于传递函数的谱表示法，其主要步骤为：（1）结合目标非高斯随机过程的概率信息确定非高斯随机过程和潜在高斯随机过程的传递函数；（2）基于目标非高斯随机过程的功率谱函数确定潜在高斯随机过程的功率谱函数；（3）结合潜在高斯随机过程的功率谱函数，采用高斯随机过程模拟方法得到相应的时程样本；（4）基于已确定的传递函数，将潜在高斯随机过程的时程样本映射得到非高斯随机过程的时程样本。传递函数的谱表示法主要包括三部分内容：传递函数的确定、潜在高斯随机过程的功率谱函数确定以及高斯随机过程模拟。根据已知非高斯随机过程的概率信息，确定传递函数的方法可分为两种：第一种是基于非高斯随机过程的概率密度函数建立无记忆非线性传递函数^[4-9]；第二种是基于非高斯随机过程的统计信息建立多项式传递函数^[10-11]。由于非高斯随机过程的统计信息可由相应的概率密度函数得到，所以第二种传递函数是第一种传递函数的近似形式。由传递函数的谱表示法的主要步骤可知，潜在高斯随机过程的功率谱函数是高斯随机过程模拟的前提条件，也是影响高斯随机过程模拟的关键因素，所以在现有的非高斯随机过程模拟研究中大多关注潜在高斯随机过程的功率谱函数计算^[4-12]，且主要分为频谱迭代更新^[4-8]和频谱修正^[9-12]。频谱迭代更新主要是经验性地调整，并不能保证所有非线性传递的收敛性^[15]，特别是对于非平稳潜在高斯随机过程的功率谱函数的计算，其计算效率和计算精度通常会受到时间点数目以及收敛准则的影响。相对而言，频谱修正通过修正不兼容的功率谱函数，从而直接获得潜在高斯随机过程的功率谱函数，是一种高效的非迭代方法。在已知高斯随机过程功率谱的前提下，谱表示法是将一系列谐和分量进行叠加，得到高斯随机过程样本。当模拟多点高斯随机过程时，通常需要对互功率谱矩阵进行Cholesky分解^[16]。当模拟随机过程的点数和频率离散数目较多时，互功率谱矩阵分解效率将会降低。针对此问题，有研究者通过插值方法^[17-20]减少互功率谱矩阵分解的次数，或通过波数-频率谱的方式^[21-24]避免互功率谱矩阵分解。在已知传递函数的前提下，潜在高斯随机过程的功率谱函数确定以及高斯随机过程模拟是非高斯随机过程模拟的重要组成部分，也是影响模拟效率的决定因素。关于多点高斯随机过程模拟中的Cholesky分解问题已发展多种高效解决方法，本文则主要关注潜在高斯随机过程的功率谱函数的高效计算方法。同时，工程结构在服役期间所遭受的荷载不仅具有非高斯概率特性，而且其概率特性具有时变性^[25-27]，所以采用非平稳非高斯随机过程对工程荷载进行描述更为合理。因此，本文将针对非平稳非高斯随机过程模拟开展研究。

值得指出的是，由于Mehler公式能够精确且高效地应用于Nataf变换中等效相关系数的计算^[28-30]，所以文献[9]通过Mehler公式构建相关函数方程，并通过求解方程确定潜在高斯随机过程的等效相关函数，虽然该方法能够有效避免求解相关函数时所遇到的数值积分问题，但对于非平稳非高斯随机过程模拟中的每一个时间点和时间延迟都需要求解方程，这必然会导致计算效率下降。为此，本文在文献[9]的基础上，通过Mehler公式建立时变相关函数方程，并针对时变相关函数方程求解所遇到的效率问题，提出一种快速插值求解方法，并以包含较少随机变量的随机过程模拟为目标，基于随机谐和函数提出非平稳非高斯随机过程快速模拟方法。

1　多点非平稳非高斯随机过程模拟

收起

模拟具有特定目标功率谱和概率密度函数的多点非平稳非高斯随机过程，采用无记忆非线性变换的思路，在时间t，将多点非平稳高斯向量随机过程Y（t）=｛Y₁（t）， Y₂（t），…， Y_n（t）｝映射得到多点非平稳非高斯向量随机过程Z（t）=｛Z₁（t），Z₂（t），…，Z_n（t）｝，即

Z j (t) = F j - 1 {Φ j [Y j (t), t], t} ≜ H [Y j (t), t]

(1)

式中，Z_j（t）∈Z（t）（j∈｛1，2，…，n｝）为指定时变边缘概率密度函数的非平稳非高斯过程；Y_j（t）∈Y（t）为零均值且方差为1的潜在非平稳高斯过程；F_j（·，t）和Φ_j（·，t）分别表示时间t对应的非高斯边缘累计分布函数和高斯边缘累计分布函数；

F j - 1

（·，t）为F_j（·，t）的逆函数；H（·，t）为转换函数。

由传递函数的谱表示法进行非高斯随机过程模拟的主要步骤可知，在式（1）建立的非高斯随机过程和潜在高斯随机过程的传递函数基础之上，精确且高效地求解潜在高斯随机过程的功率谱S^Y（ω，t）是非高斯随机过程模拟的主要步骤之一。为此，本文结合Mehler公式建立潜在高斯随机过程的功率谱S^Y（ω，t）求解方法，并针对求解S^Y（ω，t）时所遇到的计算效率问题进一步提出快速插值求解方法，进而提出多点非平稳非高斯随机过程快速模拟方法。

1.1　潜在高斯随机过程功率谱计算

基于时变Wiener-Khinchin定理^[6]，建立非平稳非高斯随机过程的时变功率谱S

j, k Z

（ω，t）（j=1，2，…，n；k=1，2，…，n）和时变自相关函数R

j, k Z

（t，τ）之间的转换关系，即

R j, k Z (t, τ) = 1 2 π ∫ - ∞ + ∞ [S j, k Z (ω, t) S j, k Z (ω, t + τ)] 12 e i ω τ d ω

(2)

式中，i²=-1；ω为频率；τ为时间延迟。

非高斯随机过程Z（t）和对应潜在高斯随机过程Y（t）的概率信息转换关系可以表示为：

ρ j, k Z (t, τ) = R j, k Z (t, τ) σ j, j Z (t) σ k, k Z (t) = E [(Z j (t) - μ Z, j (t)) (Z k (t + τ) - μ Z, k (t + τ))] σ j, j Z (t) σ k, k Z (t) = ∫ - ∞ + ∞ ∫ - ∞ + ∞ F j - 1 {Φ j [y j (t), t], t} - μ Z, j (t) σ j, j Z (t) · F k - 1 {Φ k [y k (t + τ), t], t} - μ Z, k (t + τ) σ k, k Z (t) · φ (y j (t), y k (t + τ), ρ j, k Y (t, τ)) d y j d y k

(3)

式中，ρ

j, k Z

（t，τ）为Z_j（t）和Z_k（t）的相关函数；E[∙]表示期望运算；μ_Z_，_A（t）和σ

A, A Z

（t）分别为Z_A（t）的时变均值和标准差，A=j，k；

φ (y j (t), y k (t + τ), ρ j, k Y (t, τ))

为联合高斯密度函数，可表示为：

φ (y j (t), y k (t + τ), ρ j, k Y (t, τ)) = 1 2 π 1 - (ρ j, k Y (t, τ)) 2 ⋅ e x p (- y j 2 (t) + y k 2 (t + τ) - 2 ρ j, k Y (t, τ) y j (t) y k (t + τ) 2 (1 - (ρ j, k Y (t, τ)) 2))

(4)

式中，ρ

j, k Y

（t，τ）为Y_j（t）和Y_k（t）的相关函数。

根据Mehler公式^[28-30]，式（3）可转化为关于ρ

j, k Y

（t，τ）的多项式函数，即

ρ j, k Z (t, τ) ≈ ∑ m = 1 M 1 m! I j, m (t) I k, m (t + τ) (ρ j, k Y (t, τ)) m ≜ T ρ (t, τ, ρ j, k Y (t, τ))

(5)

式中，T_ρ（t，τ，·）为在时刻t和时间延迟τ的高次多项式函数；M为Hermite多项式截断项的阶次。

I A, m (t) = ∫ - ∞ + ∞ F A - 1 [Φ (y A (t), t), t] - μ Z, A (t) σ A, A Z (t) ⋅ H m (y A) φ (y A) d y A; A = j, k

(6)

式中，H_m（·）为第m阶的Hermite多项式（m=1，2，…），具体表达式如下：

{H 1 (z i) = z i H 2 (z i) = z i 2 - 1 ⋮ H m + 1 (z i) = z i H m (z i) - m H m - 1 (z i)

(7)

当式（6）采用数值积分时，存在：

I A, m (t) = w G H, l (t) σ A, A Z π {F A - 1 [Φ (2 y G H, l (t), t), t] - μ Z, A (t)} H m (2 y G H, l (t))

(8)

式中，y_GH_，l（t）和w_GH，_l（t）为对应t时刻的高斯积分点和权重；d为数值积分的求积节点数。为保证Mehler公式求解方法的精度，d和M的取值依据分布函数不同将会有所差异，具体数值的确定参考文献[30]。

通过式（5）求解M次代数方程即可确定等效相关函数ρ

j, k Y

（t，τ），随后潜在高斯随机过程的功率谱函数S

j, k Y

（ω，t）由下式计算：

S j, k Y (ω, t) = ∫ - ∞ + ∞ [ρ j, k Y (t, τ)] e - i ω τ d τ

(9)

由式（5）可知，为确定潜在高斯随机过程的相关函数ρ

j, k Y

（t，τ），需求解高次多项式方程。显然，求解高次多项式方程的次数与时间点和时间延迟数目有关，这必然会影响模拟方法的计算效率。为提高计算效率，将在下文提出潜在高斯随机过程功率谱快速计算方法。

1.2　潜在高斯随机过程功率谱快速计算方法

对于指定的非高斯随机过程，高次多项式函数T_ρ（t，τ，·）可通过I_A_，_m（t）确定，则称通过ρ

j, k Z

（t，τ）求解相应ρ

j, k Y

（t，τ）为正问题，对应的逆问题通过ρ

j, k Y

（t，τ）计算ρ

j, k Z

（t，τ）。显然，计算多项式函数的耗时远低于求解相应次数的非线性方程。为此，本节在求解逆问题的基础上，结合插值方法快速计算ρ

j, k Y

（t，τ）。

由于等效相关函数ρ

j, k Y

（t，τ）的取值范围在区间[-1，1]，所以采用间隔2/（s-1）将区间均分为s个确定的计算点

{ρ ¯ Y 1, ρ ¯ Y 2, ⋅ ⋅ ⋅, ρ ¯ s − 1 Y, ρ ¯ s Y} ≜ ρ ¯ Y

，即

ρ ̲ l Y = - 1 + (l - 1) 2 s - 1; l = 1,2, …, s

(10)

将

ϱ ¯ Y

中的数值依次代入由I_A_，_m（t）确定的多项式函数T_ρ（t，τ，·）中，可确定相应数值

ρ ¯ Z ≜ {ρ ¯ Z 1, ρ ¯ 2 Z, ⋯, ρ ¯ s − 1 Z, ρ ¯ s Z}

。由数值

ρ ¯ Y

与

ϱ ¯ Z

可确定在时刻t和时间延迟τ时

ρ j, k Y (t, τ)

与

ρ j, k Z (t, τ)

的对应关系（如图1中的蓝色线）。在本文中，s=1000。

通过已确定的非高斯随机过程和潜在高斯随机过程的相关函数对应关系，采用插值的方式计算得到ρ

j, k Z

（t，τ）所对应的ρ

j, k Y

（t，τ）（如图2所示），其中插值函数可采用三次样条插值，即MATLAB软件中的“interp1”命令。由于ρ

j, k Z

（t，τ）与ρ

j, k Y

（t，τ）在固定时间t和时间延迟τ的情况下保持单调递增关系^[31]，所以ρ

j, k Z

（t，τ）和ρ

j, k Y

（t，τ）保持一一对应关系。

当潜在高斯随机过程的等效相关函数确定后，采用式（9）可得到对应演变功率谱。对于无记忆非线性平移存在的概率信息可能出现不相容问题以及非正定问题，潜在高斯随机过程的演变功率谱需进行谱修正，具体的修正方式参考文献[9]。

1.3　基于随机谐和函数的多点非平稳高斯随机过程模拟

当潜在非平稳高斯随机过程Y（t）的功率谱函数确定后，可采用非平稳高斯随机过程模拟方法生成时程样本。当随机过程应用于工程结构的随机响应分析时，通常需要以较少的随机变量生成时程样本，例如概率密度演变法。孙伟玲等^[32]、CHEN等^[33-35]采用少量的随机谐和分量对高斯随机过程进行模拟且能够较好地拟合目标时变功率谱。因此，本文采用随机谐和函数进行多点非平稳高斯随机过程模拟。

根据1.2节中确定的潜在非平稳高斯随机过程的功率谱函数S

j, k Y

（ω，t），Y（t）的功率谱矩阵S^Y（ω，t）可由S

j, k Y

（ω，t）表示为：

S Y (ω, t) = [S 1,1 Y (ω, t) S 1,2 Y (ω, t) ⋯ S 1, n Y (ω, t) S 2,1 Y (ω, t) S 2,2 Y (ω, t) ⋯ S 2, n Y (ω, t) ⋮ ⋮ ⋮ S n, 1 Y (ω, t) S n, 2 Y (ω, t) ⋯ S n, n Y (ω, t)]

(11)

利用Cholesky分解，可将S^Y（ω，t）分解为：

S Y (ω, t) = H (ω, t) [H * (ω, t)] T

(12)

式中，上标“T”为转置符号；H^*（ω，t）为H（ω，t）的复共轭矩阵，H（ω，t）可表示为：

H (ω, t) = [H 1,1 (ω, t) 0 ⋯ 0 H 2,1 (ω, t) H 2,2 (ω, t) ⋯ 0 ⋮ ⋮ ⋮ H n, 1 (ω, t) H n, 2 (ω, t) ⋯ H n, n (ω, t)]

(13)

基于随机谐和函数，非平稳高斯随机过程Y_j（t）表示为：

Y j (t) = ∑ m = 1 j ∑ l = 1 N | A j, m (ω m l, t) | c o s (ω m l t + ϕ m l)

(14)

式中，N为谐和分量数目，可由文献[34]确定；ϕ_ml为[0，2π]区间的独立同分布的均匀变量，其密度函数相应为：

p ϕ m l = 1 2 π

(15)

| A j, m (ω m l, t) | = | H j, m (ω m l, t) | ω m l (p) - ω m (l - 1) (p) π

(16)

随机频率ω_ml的密度函数

p ω m l (ω)

为：

p ω m l (ω) = {1 ω m l (p) - ω m (l - 1) (p), ω ∈ [ω m (l - 1) (p), ω m l (p)] 0, ω ∉ [ω m (l - 1) (p), ω m l (p)]

(17)

式中，ω

m l (p)

为[ω_m_L，ω_m_U]内的任意内插点，其中ω_m_L和ω_m_U分别为随机过程Y_l（t）的截止频率下限和上限。

由式（14）获得潜在高斯随机过程后，结合式（1），即可获得非高斯随机过程。

1.4　多点非平稳非高斯随机过程模拟步骤

关于本文提出的非平稳非高斯随机过程模拟方法的计算步骤总结如下：

（1）根据非高斯随机过程Z（t）的功率谱S

j, k Z

（ω，t）和式（3），确定非高斯随机过程的相关函数ρ

j, k Z

（t，τ）；

（2）结合Z（t）的时变概率信息，确定不同时刻t处d和M的取值，并通过

ρ ¯ Y

与

ϱ ¯ Z

的离散点建立ρ

j, k Z

（t，∙）和ρ

j, k Y

（t，∙）的映射关系；

（3）基于已建立的ρ

j, k Z

（t，∙）和ρ

j, k Y

（t，∙）的映射关系，通过插值方法计算ρ

j, k Z

（t，τ）所对应的ρ

j, k Y

（t，τ），进而计算潜在非平稳高斯随机过程的功率谱函数S

j, k Y

（ω，t）；

（4）根据潜在高斯随机过程功率谱函数S

j, k Y

（ω，t）和式（14）进行非平稳高斯随机过程模拟，得到相应时程样本过程；

（5）通过式（1）和非平稳高斯时程样本过程，得到具有特定目标功率谱和概率信息的非平稳非高斯时程样本过程。

2　算例验证

收起

为验证本文所提非平稳非高斯随机过程模拟方法（记为建议法）的有效性，本节将通过两个算例对其进行验证，其中算例1为单点均匀调制非高斯随机过程模拟，算例2为多点非均匀调制非高斯随机过程模拟。为对比分析，本文同时考察文献[6]提出的潜在高斯随机过程的功率谱密度函数迭代求解方法，并将其与谱表示法^[16]和无记忆非线性变换^[4]结合，获得相应非高斯随机过程（记为迭代法）。

随机过程Z（t）的时变均值μ_Z（t）、时变二阶矩m_Z_，2（t）、时变三阶矩m_Z_，3（t）、时变四阶矩m_Z_，4（t）、时变偏度sw_Z（t）和时变峰度ks_Z（t）的计算表达式如附录所示。

本文所有程序均在处理器为i5-12600 K、主频为3.70 GHz和内存为32 GB的计算机上运行。

2.1　单点均匀调制非高斯随机过程模拟

对于单点均匀调制非高斯随机过程Z（t），其边缘密度函数为Gamma分布，且均值μ_Z（t）、二阶矩m_Z，2（t）、偏度sw_Z（t）和四阶矩m_Z，4（t）的目标值如图3所示。Z（t）可由下式表示：

Z (t) = g (t) X (t)

(18)

式中，X（t）为平稳非高斯随机过程；g（t）为均匀调制函数，表示为：

g (t) = e - 0.25 t - e - 0.5 t 0.25

(19)

Z（t）的演变功率谱密度函数S^Z（ω，t）为：

S Z (ω, t) = g 2 (t) S X (ω)

(20)

式中，S^X（ω）为X（t）的功率谱密度函数，表示为：

S X (ω) = ω g 4 + 4 ξ g 2 ω g 2 ω 2 (ω 2 - ω g 2) 2 + 4 ξ g 2 ω g 2 ω 2 S 0

(21)

式中，ω_g=12.15 rad/s，ξ_g=0.65，S₀=0.07 m²∙s^-3。计算中，截止频率的上、下限分别为ω_U =100 rad/s和ω_L=0，频率离散数目N_w为1000，时间延迟离散数目N_τ为256，时间离散数目N_t为1024。

对于均匀调制非高斯随机过程模拟，通常可对平稳非高斯随机过程X（t）进行模拟得到样本随机过程，再结合式（18），即可得到Z（t）的样本随机过程，其中结合文献[30]，d和M的取值分别为11和3。为验证Mehler公式求解相关函数的精度，图3分别绘制了由式（3）通过积分求解（记为标准解）、Mehler公式求解以及插值求解Y（t）的相关函数。由图4可知，通过Mehler公式求解得到Y（t）的相关函数能够与标准解在整个时间延迟上吻合，这表明Mehler公式求解具有较高的计算精度，同时，插值求解的结果与Mehler公式求解的结果完全吻合，这表明插值求解能够精确确定潜在高斯随机过程的相关函数。随后，基于1.1节和1.2节，分别通过Mehler公式求解和插值求解，得到X（t）对应的潜在高斯随机过程Y（t）的功率谱密度函数S^Y（ω），如图4所示。同时，由迭代求解得到的潜在高斯随机过程的功率谱密度函数绘制于图4。由图4可知，通过插值方法和Mehler公式求解得到的S^Y（ω）能够与标准解吻合，而通过迭代求解得到的S^Y（ω）虽然能够与标准解相近，但在峰值处却表现出明显差异。为进一步考察Mehler公式求解、插值求解和迭代求解的计算精度，通过下式计算上述三种方法得到的S^Y（ω）与目标解的差异ε_G，即

ε G = ∫ 0 + ∞ [S T Y (ω) - S E Y (ω)] 2 d ω ∫ 0 + ∞ [S T Y (ω)] 2 d ω × 100 %

(22)

式中，S

E Y

（ω）为Mehler公式求解、插值求解和迭代求解得到的潜在高斯随机过程的功率谱，S

T Y

（ω）为相应的标准解。

由式（22）得到Mehler公式求解、插值求解和迭代求解与目标解的差异分别为0.0943%、0.0943%和9.97%。显然，Mehler公式求解和插值求解能够更为精确地计算潜在高斯随机过程的功率谱。需指出的是，通过Mehler公式求解S^Y（ω）的时间为68.25 s，通过插值求解S^Y（ω）的时间为0.51 s，而迭代求解花费9.06 s计算S^Y（ω）。其中在插值求解中，0.19 s用于确定

ρ ¯ Y

与

ϱ ¯ Z

的关系，0.32 s用于拟合和计算S^Y（ω）。同时，依据时间复杂度的概念^[36]，插值求解需要在每一个时间延迟点处计算潜在高斯随机过程的相关函数，共需循环计算N_τ次，随后需要循环N_w次计算潜在高斯随机过程的功率谱函数，则插值求解的时间复杂度为O（N_τ+N_w）。迭代求解在每一次迭代循环中采用N_τ次计算迭代非高斯随机过程的相关函数，N_w次计算迭代非高斯随机过程的功率谱函数。在本算例中，共采用10次循环判断迭代非高斯随机过程的功率谱函数收敛目标非高斯随机过程的功率谱函数，则迭代求解的时间复杂度为O[10∙（N_τ+N_w）]。显然，迭代求解的时间复杂度明显高于插值求解的时间复杂度，进一步表明，插值求解算法程序的高效性。上述结果表明插值求解能够高效且精确地计算潜在高斯随机过程的功率谱密度函数S^Y（ω）。

在获得S^Y（ω）后，采用14个谐和分量进行高斯随机过程模拟，得到1000条Y（t）的样本过程，再通过无记忆非线性变换得到对应X（t）的时程样本，最终通过式（18）获得Z（t）的样本过程。由附录分别计算Z（t）的时变统计矩，并绘制于图5。由图5可知，由建议法得到非高斯随机过程的前四阶矩能够与目标值吻合，而迭代法得到的统计矩却与目标值有差异。图6绘制了3条非平稳非高斯随机过程Z（t）的时程样本。从图6中可以看出，每条时程样本在整个时域上表现出明显的非平稳特性。图7分别绘制建议法和迭代法在t =3 s时的概率密度函数。由图7可知，建议法的密度函数能够和目标密度函数在整个区间吻合，而迭代法的密度函数在峰点处与目标密度函数有差异。图8分别绘制了由建议法和迭代法得到的时变功率谱。显然，建议法得到非高斯随机过程的时变功率谱与目标谱在整个频域和时域上具有一致性，迭代法得到的时变功率谱幅值却与目标谱有一定差异。需指出的是，建议法和迭代法生成一条Z（t）的样本过程分别耗时0.513和9.17 s。由于建议法和迭代法在高斯随机过程模拟时采用矩阵运算且两种方法循环次数一致，所以本算例没有对比两种方法在高斯随机过程模拟时的时间复杂度。本算例的计算结果表明，建议法相对迭代法能够有效地模拟具有特定目标功率谱的非高斯随机过程。

2.2　多点非均匀调制非高斯随机过程模拟

边缘密度函数为对数正态分布的非均匀调制非高斯向量随机过程Z（t）=｛Z₁（t），Z₂（t），…，Z₁₀（t）｝，其均值依次为3、5、7、9、11、13、15、17、19、21，时变二阶矩m_Zl_，2（t）（l = 1，2，…，10）如图9所示，Z（t）的其他分布参数和统计信息可通过相应的密度函数计算得到，且部分目标时变统计矩m_Z2，2（t）、m_Z2，3（t）和ks_Z2（t）如图10所示。Z_i（t）和Z_j（t）的互功率谱S

i, j Z

（ω，t）为：

S i j Z (ω, t) = A i (ω, t) A j (ω, t) S i j X (ω); i = 1,2, …, 10; j = 1,2, …, 10; i ≠ j

(23)

式中，A_l（ω，t）（l = 1，2，…，10）为调制函数并采用相同的形式，即

A l (ω, t) = e x p (- a t) - e x p (- (c ω + b) t) e x p (- a t d) - e x p (- (c ω + b) t d)

(24)

式中，a=0.25 s^-1，b=0.251 s^-1，c=0.05。

t d = l n (c ω + b) - l n a c ω + (b - a)

(25)

S i j X (ω) = S i i X (ω) S j j X (ω) Ψ i j X (ω)

(26)

其中：

S j j X (ω) = ω g 4 + 4 ξ g 2 ω g 2 ω 2 (ω 2 - ω g 2) 2 + 4 ξ g 2 ω g 2 ω 2 ⋅ ω 4 (ω 2 - ω f 2) + 4 ξ f 2 ω f 2 ω 2 ⋅ C j

(27)

Ψ i j X (ω) = e x p [- a (ω) d i, j]

(28)

{a (ω) = 0.1678 × 10 - 4 ω 2 + 0.1219 × 10 - 2 b (ω) = - 0.0055 ω + 0.7674

(29)

式中，常数ω_g、ω_f、ξ_g、ξ_f分别取为8.5 rad/s、0.85 rad/s、0.75、0.75；C_j（j=1，2，…，10）依次为0.08、0.10、0.12、0.14、0.16、0.18、0.20、0.22、0.24、0.26 m²/s³；d_i，j计算如下：

d i, j = 300 j + (i - 1) × 100

(30)

在随机过程模拟中，截止频率的上、下限分别为ω_U=100 rad/s和ω_L=0，频率离散数目N_w为1000，时间延迟离散数目N_τ为256，时间离散数目N_t为256。

由于对数正态分布的相关系数和对应高斯分布的相关系数有解析解^[31]，所以本算例将通过解析解得到的潜在高斯随机过程的功率谱作为标准解。

基于文献[30]，d的取值恒为11，不同时刻处M的取值如图11所示。由图11可知，除Z₁（t）在特定时间段内M的取值为3，其余时刻不同随机过程所对应的M取值恒为2。根据非平稳非高斯随机过程的时变功率谱，分别由Mehler公式求解、插值求解和迭代求解得到潜在高斯随机过程的时变功率谱密度函数，以S

2,6 Y

（ω，t）为例绘制图12。由图12可知，Mehler公式求解和插值求解在整个时频域内与标准解吻合，迭代求解在低频处能够与标准解吻合，而在高频处却与标准解表现出明显差异。同时，Mehler公式求解、插值求解、迭代求解与标准解的差异可由下式计算：

ε G = ∫ 0 ∞ ∫ 0 + ∞ [S T Y (ω, t) - S E Y (ω, t)] 2 d ω d t ∫ 0 + ∞ [S T Y (ω, t)] 2 d ω d t × 100 %

(31)

式中，S

E Y

（ω，t）为Mehler公式求解、插值求解、迭代求解计算得到的潜在高斯随机过程的功率谱；S

T Y

（ω，t）为潜在高斯随机过程的功率谱标准解。

由式（31）得到Mehler公式求解、插值求解和迭代求解与目标解的差异分别为0.7306%、0.7306%和32.7158%。Mehler公式求解、插值求解和迭代求解得到S

2,6 Y

（ω，t）所花费的时间依次为77824、2.15和2821.65 s。同时，插值求解在每一个时间点处需进行N_τ次潜在高斯随机过程的相关函数计算和N_w次潜在高斯随机过程的功率谱函数计算，迭代求解程序在每一次迭代循环中的每一个时间点处涉及N_τ次迭代非高斯随机过程的相关函数计算和N_w次迭代非高斯随机过程的功率谱函数计算，本算例共采用15次迭代循环，则插值求解和迭代求解的时间复杂度依次为O[N_t∙（N_τ+N_w）]和O[15∙N_t∙（N_τ+N_w）]。显然，插值求解能够高效且精确地计算潜在高斯随机过程的功率谱函数。

在确定潜在高斯随机过程的时变功率谱后，采用76个谐和分量进行非平稳高斯随机过程模拟方法，得到潜在高斯随机过程的1000条时程样本，并结合式（1），即可获得非高斯随机过程的时程样本。图13分别绘制了Z₂（t）、Z₅（t）和Z₉（t）的时程样本。从图13中可以明显看出，时程样本在整个时间域上表现出明显的非平稳特性；同时，不同点处的时程样本幅值有较大差异。基于建议法和迭代法得到的样本过程，由附录分别计算时变统计矩M_Z_2，2（t）、M_Z_2，3（t）和ks_Z₂（t），并绘制于图10。由图10可知，建议法和迭代法的时变统计矩能够与目标值在整个时间域上吻合，且具有明显的非平稳特性。图14给出了Z₂（t）在t=8.6 s时的概率密度函数。由图14可知，建议法和迭代法所对应的概率密度函数能够和目标密度函数匹配。通过小波方法^[37]，估计非平稳非高斯样本随机过程的时变功率谱，其中图15分别给出S

2,2 Z

（ω，t）的目标功率谱和估计得到的时变功率谱。由图15可知，建议法和迭代法所对应的时变功率谱在时域和频域上能够与目标时变功率谱整体保持一致。在时变功率谱的初始时刻，目标功率谱和估计功率谱存在一定差异，该差异的主要原因是小波估计的局限性所导致^[34]。建议法和迭代法生成一条Z₂（t）的样本过程分别耗时52.22和91981.37 s。显然，建议法相对迭代法能够有效地模拟具有特定目标功率谱的非高斯随机过程。

为进一步验证建议法模拟非高斯随机过程的高效性，通过对非高斯向量随机过程的维数nn依次取为2、3、4、5、6、7、8、9、10，并计算生成一条随机过程的时间，如图16所示。由图16可知，建议法的耗时随着维数nn的增加而增长，但整体耗时仍处于可接受的范围。

3　结论

收起

本文针对多点非平稳非高斯随机过程，采用插值求解和随机谐和函数，提出一种快速模拟方法，结论如下：

（1）相对于Mehler公式求解和迭代求解潜在高斯随机过程的功率谱，插值求解在保持计算精度的同时能够提高计算效率，且在计算多点潜在高斯随机过程的时变功率谱时，插值求解的优势更为明显。

（2）通过随机谐和函数，能够以较少的谐和分量进行非平稳高斯随机过程模拟。

（3）通过非平稳非高斯随机过程快速模拟方法得到样本过程，能够与目标功率谱、目标统计矩和目标密度函数吻合。

（4）对于非高斯随机过程模拟，本文所提插值求解方法不仅可以与随机谐和函数相结合，还可以与现有高效地高斯随机过程模拟相结合。

（5）在多点非平稳非高斯随机过程模拟中，如何高效地分解功率谱矩阵是未来研究的重点。

基金

收起

河南省高等学校重点科研项目(24A560007)
国家自然科学基金资助项目(52102096)
国家自然科学基金资助项目(51678092)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

杨庆山，田玉基. 地震地面运动及其人工合成[M]. 北京：科学出版社，2014:20-50.

YANG

Qingshan

，TIAN

Yuji

. Earthquake Ground Motions & Artificial Generation[M]. Beijing：Science Press，2014:20-50.

[2]

黄国庆，彭留留，廖海黎，等.普立特大桥桥位处山区风特性实测研究[J]. 西南交通大学学报，2016,51(2):349-356.

HUANG

Guoqing

，PENG

Liuliu

，LIAO

Haili

，et al.Field measurement study on wind characteristics at Puli Great Bridge site in mountainous area[J]. Journal of Southwest Jiaotong University，2016,51(2):349-356.

[3]

吴家驹，张鹏飞，胡亚冰. 非高斯随机振动的分析基础[J]. 强度与环境，2018,45(2):1-8.

Jiaju

，ZHANG

Pengfei

，HU

Yabing

. Analytical basis for the synthesis of non-Gaussian random vibration[J]. Structure & Environment Engineering，2018,45(2):1-8.

[4]

GRIGORIU

. Applied Non-Gaussian Processes[M]. New Jersey：Prentice Hall，1995.

[5]

FERRANTE

F J

，ARWADE

S R

，GRAHAM-BRADY

L L

. A translation model for non-stationary，non-Gaussian random processes[J]. Probabilistic Engineering Mechanics，2005,20(3):215-228.

[6]

SHIELDS

M D

，DEODATIS

. Estimation of evolutionary spectra for simulation of non-stationary and non-Gaussian stochastic processes[J]. Computers & Structures，2013,126:149-163.

[7]

Y X

，GAO

Y F

，ZHANG

，et al.Simulation of spatially varying non-Gaussian and nonstationary seismic ground motions by the spectral representation method[J]. Journal of Engineering Mechanics，2018,144(1):04017143.

[8]

，PENG

L L

，JIANG

，et al.Simulation of stationary non-Gaussian stochastic vector processes using an eigenvalue-based iterative translation approximation method[J]. Mechanical Systems and Signal Processing，2022,175:109128.

[9]

Z L

，LIU

R Y

，FAN

W L

，et al.Simulation of stationary Gaussian/non-Gaussian stochastic processes based on stochastic harmonic functions[J]. Probabilistic Engineering Mechanics，2021,66:103141.

[10]

FAN

W L

，TIAN

，SHENG

X Q

. Simulation of a non-Gaussian stochastic process based on a combined distribution of the UHPM and the GBD[J]. Probabilistic Engineering Mechanics，2023,72:103438.

[11]

Z H

，ZHAO

，ZHANG

X Y

，et al.Simulating stationary non-Gaussian processes based on unified Hermite polynomial model[J]. Journal of Engineering Mechanics，2020,146 (7):04020067.

[12]

李锦华，李春祥，蒋磊，等.多变量非高斯脉动风压模拟的非迭代法[J].上海交通大学学报，2017,51(4):462-469.

Jinhua

，LI

Chunxiang

，JIANG

Lei

，et al.A non-iterative approach for simulating multivariate non-Gaussian fluctuating wind pressure[J]. Journal of Shanghai Jiao Tong University，2017,51(4):462-469.

[13]

PHOON

K K

，HUANG

H W

，QUEK

S T

. Simulation of strongly non-Gaussian processes using Karhunen-Loeve expansion[J]. Probabilistic Engineering Mechanics，2005,20(2):188-198.

[14]

SHIELDS

M D

，KIM

. Simulation of higher-order stochastic processes by spectral representation[J]. Probabilistic Engineering Mechanics，2017,47:1-15.

[15]

PUIG

，AKIAN

J L

. Non-Gaussian simulation using Hermite polynomials expansion and maximum entropy principle[J]. Probabilistic Engineering Mechanics，2004,19(4):293-305.

[16]

SHINOZUKA

，JAN

C M

. Digital simulation of random processes and its applications[J]. Journal of Sound and Vibration，1972,25(1):111-128.

[17]

ZHAO

，JIANG

，PENG

L L

，et al.Fast simulation of nonstationary wind velocity fields by proper orthogonal decomposition interpolation[J]. Journal of Wind Engineering and Industrial Aerodynamics，2021,219:104798.

[18]

TAO

T Y

，WANG

，ZHAO

K Y

. Efficient simulation of fully non-stationary random wind field based on reduced 2D Hermite interpolation[J]. Mechanical Systems and Signal Processing，2021,150:107265.

[19]

杨雄骏，黄金山，张建国，等.基于插值与降维方法的输电塔线体系随机脉动风场有效模拟[J]. 振动与冲击，2021,40(9):77-83.

YANG

Xiongjun

，HUANG

Jinshan

，ZHANG

Jianguo

，et al.Effective simulation of stochastic fluctuating wind field of transmission tower-line system based on interpolation and dimension reduction method[J]. Journal of Vibration and Shock，2021,40(9):77-83.

[20]

JIANG

，ZHAO

，PENG

L L

，et al.Fast simulation of fully non-stationary wind fields using a new matrix factorization assisted interpolation method[J]. Mechanical Systems and Signal Processing，2022,172:108973.

[21]

宋玉鹏，陈建兵，彭勇波. 基于波数-频率联合演变功率谱的一维空间非均匀脉动风场模拟[J]. 工程力学，2019,36(2):205-214.

SONG

Yupeng

，CHEN

Jianbing

，PENG

Yongbo

. Simulation of nonhomogeneous fluctuating wind field in one-dimensional space by evolutionary wavenumber-frequency joint power spectrum[J]. Engineering Mechanics，2019,36(2):205-214.

[22]

YANG

X J

，LEI

. Efficient simulation of wind fields based on the factorization of wavenumber-frequency joint spectrum[J]. Computer-Aided Civil and Infrastructure Engineering，2022,37(3):370-385.

[23]

ZHAO

，LU

Z H

，ZHAO

Y G

. Simulating multivariate stationary non-Gaussian process based on wavenumber-frequency spectrum and unified Hermite polynomial model[J]. Probabilistic Engineering Mechanics，2022,69:103272.

[24]

，XU

. A PDF discretization scheme in wavenumber-frequency joint spectrum for simulating multivariate random fluctuating wind fields[J]. Probabilistic Engineering Mechanics，2023,72:103422.

[25]

GUPTA

I D

，TRIFUNAC

M D

. A note on the nonstationarity of seismic response of structures[J]. Engineering Structures，2000,22(11):1567-1577.

[26]

CHEN

X Z

. Analysis of alongwind tall building response to transient nonstationary winds[J]. Journal of Structural Engineering，2008,134(5):782-791.

[27]

FALSONE

，SETTINERI

. A method for the random analysis of linear systems subjected to non-stationary multi-correlated loads[J]. Probabilistic Engineering Mechanics，2011,26(3):447-453.

[28]

KIBBLE

W F

. An extension of a theorem of Mehler’s on Hermite polynomials[J]. Mathematical Proceedings of the Cambridge Philosophical Society，1945,41(1):12-15.

[29]

范文亮，杨朋超，李正良. 基于Mehler公式的等效相关系数求解技术[J]. 同济大学学报（自然科学版），2016,44(6):838-844.

FAN

Wenliang

，YANG

Pengchao

，LI

Zhengliang

. A technique for solution of equivalent correlation coefficients based on the Mehler’s formula[J]. Journal of Tongji University （Natural Science），2016,44(6):838-844.

[30]

Y F

，FAN

W L

，ZHANG

J Y

，et al.Investigation of equivalent correlation coefficient based on the Mehler’s formula[J]. Engineering Computations，2019,36(4):1169-1200.

[31]

LIU

P L

，DER KIUREGHIAN

. Multivariate distribution models with prescribed marginals and covariances[J]. Probabilistic Engineering Mechanics，1986,1(2):105-112.

[32]

孙伟玲，陈建兵，李杰. 随机过程的第二类随机谐和函数表达[J]. 同济大学学报（自然科学版），2011,39(10):1413-1419.

SUN

Weiling

，CHEN

Jianbing

，LI

Jie

. Stochastic harmonic functions of second kind for spectral representations[J]. Journal of Tongji University （Natural Science），2011,39(10):1413-1419.

[33]

CHEN

J B

，SUN

W L

，LI

，et al.Stochastic harmonic function representation of stochastic processes[J]. Journal of Applied Mechanics，2012,80(1):011001.

[34]

CHEN

J B

，KONG

，PENG

Y B

. A stochastic harmonic function representation for non-stationary stochastic processes[J]. Mechanical Systems and Signal Processing，2017,96:31-44.

[35]

CHEN

J B

，COMERFORD

，PENG

Y B

，et al.Reduction of random variables in the stochastic harmonic function representation via spectrum-relative dependent random frequencies[J]. Mechanical Systems and Signal Processing，2020,141:106718.

[36]

王红梅，皮德常.数据结构：从概念到C实现[M]. 北京：清华大学出版社，2017.

WANG

Hongmei

，PI

Dechang

. Data Structure：From Concept to C Implementation[M]. Beijing：Tsinghua University Press，2017.

[37]

KONG

，SPANOS

P D

，LI

，et al.Response evolutionary power spectrum determination of chain-like MDOF non-linear structural systems via harmonic wavelets[J]. International Journal of Non-Linear Mechanics，2014,66:3-17.

2025年第38卷第8期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.16385/j.cnki.issn.1004-4523.202309035

接收时间：2023-09-13
首发时间：2026-02-09
出版时间：2025-08-10

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2023-09-13
修回日期：2023-11-27

基金

河南省高等学校重点科研项目(24A560007)

国家自然科学基金资助项目(52102096)

国家自然科学基金资助项目(51678092)

作者信息

^1.河南科技大学工程力学系，河南洛阳 471000

^2.重庆大学建筑力学系，重庆 400045

^3.西安机电信息研究所，陕西西安 710065

通讯作者:

虞跨海（1982—），男，博士，教授。E-mail: yukuahai@163.com

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/zdgcxb/CN/10.16385/j.cnki.issn.1004-4523.202309035

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT