Page 189 - 《振动工程学报》2025年第8期

第 8 期 盛向前，等： 非平稳非高斯随机过程插值模拟方法 1829

1. 1 潜在高斯随机过程功率谱计算 2，…），具体表达式如下：

ì H 1 ( z i )= z i
ï
基于时变 Wiener‑Khinchin 定理 ，建立非平稳 ï ï
［6］
2
ï ï ï H 2 ( z i )= z i - 1
Z
非 高 斯 随 机 过 程 的 时 变 功 率 谱 S j，k （ω，t）（j=1， í ⋮ （7）
ï ï
Z
ï
2，…，n；k=1，2，…，n）和时变自相关函数 R j，k （t，τ） ï
î
ï ï H m + 1 ( z i )= z i H m ( z i )- mH m - 1 ( z i )
之间的转换关系，即
当式（6）采用数值积分时，存在：
1 +∞ Z Z 1 iωτ
Z
2
R j，k(t，τ) = ∫ [ S j，k ( ω，t ) S j，k ( ω，t + τ )] e dω w GH，l( ) t { F A ( ) ù
-1é
2π -∞ I A，m(t) = ë Φ 2 y GH，l( ) t ，t ，t -
û
Z
（2） σ A，A π
)
式中，i =-1；ω 为频率；τ 为时间延迟。 μ Z，A( ) t } H m( 2 y GH，l( ) t （8）
2
非高斯随机过程 Z（t）和对应潜在高斯随机过
式中，y GH，l （t）和 w GH，l （t）为对应 t 时刻的高斯积分点
程 Y（t）的概率信息转换关系可以表示为：
R j，k( ) 和 权 重 ；d 为 数 值 积 分 的 求 积 节 点 数 。 为 保 证
Z
t，τ
Z
ρ j，k ( t，τ )= = Mehler 公式求解方法的精度，d 和 M 的取值依据分
σ j，j( ) t σ k，k( ) t
Z
Z
)
E[( Z j ( t )- μ Z，j( ) t ( Z k ( t + τ )- μ Z，k(t + τ ] ) ) 布函数不同将会有所差异，具体数值的确定参考文
= 献［30］。
σ j，j( ) t σ k，k( ) t 通过式（5）求解 M 次代数方程即可确定等效相
Z
Z
+∞ F j {Φ j[ y j( ) t ，t Y
-1
+∞ ]，t }- μ Z，j( ) t 关函数 ρ j，k （t，τ），随后潜在高斯随机过程的功率谱
∫ ∫
Z · Y
-∞ -∞ σ j，j( ) t 函数 S j，k （ω，t）由下式计算：
F k {Φ k[ y k(t + τ )，t ]，t }- μ Z，k(t + τ ) S j，k ( ω，t )= ∫ +∞ [ ρ j，k(t，τ)] e -iωτ dτ （9）
Y
Y
-1
Z · -∞
σ k，k( ) t
由式（5）可知，为确定潜在高斯随机过程的相关
Y （3） Y
φ( y j(t)，y k(t + τ)，ρ j，k(t，τ)) dy j dy k 函数 ρ j，k （t，τ），需求解高次多项式方程。显然，求解
Z
式中，ρ j，k （t，τ）为 Z j （t）和 Z k （t）的相关函数；E［∙］表示 高次多项式方程的次数与时间点和时间延迟数目有
Z
期望运算；μ Z，A （t）和 σ A，A （t）分别为 Z A （t）的时变均值 关，这必然会影响模拟方法的计算效率。为提高计
Y
和 标 准 差 ，A=j，k；φ( y j(t)，y k(t + τ)，ρ j，k(t，τ)) 为 算效率，将在下文提出潜在高斯随机过程功率谱快
联合高斯密度函数，可表示为： 速计算方法。
Y
φ( y j(t)，y k(t + τ)，ρ j，k(t，τ)) = 1. 2 潜在高斯随机过程功率谱快速计算方法
1
⋅ 对于指定的非高斯随机过程，高次多项式函数
1 -( ρ j，k( ) ) 2
Y
2π t，τ T ρ （t，τ，·）可通过 I A，m （t）确定，则称通过 ρ j，k （t，τ）求解
Z
æ ö 相应 ρ j，k （t，τ）为正问题，对应的逆问题通过 ρ j，k （t，τ）
Y
Y
)
)
t，τ
Y
ç ç y j ( ) t + y k(t + τ - 2ρ j，k( ) y j( ) t y k(t + τ ÷÷
2
2
ç
(
Z
exp - 2 1 -( ρ j，k( ) 2 ÷ ÷ ÷ ÷ 计算 ρ j，k （t，τ）。显然，计算多项式函数的耗时远低
ç
ç ç
Y
è t，τ ) ) ø 于求解相应次数的非线性方程。为此，本节在求解
Y
（4） 逆问题的基础上，结合插值方法快速计算 ρ j，k （t，τ）。
Y
Y
式中，ρ j，k （t，τ）为 Y j （t）和 Y k （t）的相关函数。 由于等效相关函数 ρ j，k （t，τ）的取值范围在区间
根 据 Mehler 公 式 ［28‑30］ ，式（3）可 转 化 为 关 于 ［-1，1］，所以采用间隔 2/（s-1）将区间均分为 s 个
Y
Y
Y
Y
ρ j，k （t，τ）的多项式函数，即 确定的计算点｛ρ Y1， ρ Y2，…， ρ s‒1， ρ s｝≜ρ ，即
M Y 2
Z 1 ) Y m ρ =-1 +(l - 1) ；l = 1，2，…，s （10）
ρ j，k(t，τ) ≈ ∑ I j，m(t) I k，m(t + τ ( ρ j，k(t，τ)) ≜ - l s - 1
m = 1 m！ Y
T ρ( t，τ，ρ j，k ( t，τ )) （5） 将 ρ 中的数值依次代入由 I A，m （t）确定的多项式
Y
Z
Z
函数 T ρ （t，τ，·）中，可确定相应数值 ρ ≜｛ρ Z1，ρ 2 ，…，
式中，T ρ （t， τ， ·）为在时刻 t 和时间延迟 τ 的高次多
Z
Y
Z
Z
ρ s - 1，ρ s｝。由数值 ρ 与 ρ 可确定在时刻 t 和时间延
项式函数；M 为 Hermite 多项式截断项的阶次。 迟 τ 时 ρ j，k （t，τ）与 ρ j，k （t，τ）的对应关系（如图 1 中的蓝
Y
Z
)
]
+∞ F A [ Φ( y A( ) t ，t ，t - μ Z，A( ) t
-1
∫ 色线）。在本文中，s=1000。
I A，m(t) = Z ⋅
-∞ σ A，A( ) t
通过已确定的非高斯随机过程和潜在高斯随机
H m( y A) φ( y A) dy A； A = j，k （6） 过程的相关函数对应关系，采用插值的方式计算得
Y
Z
式中 ，H m （·）为 第 m 阶 的 Hermite 多 项 式（m=1， 到 ρ j，k （t，τ）所对应的 ρ j，k （t，τ）（如图 2 所示），其中插