CN104516771A - 一种非平稳随机过程高效模拟方法 - Google Patents

Abstract

本发明属于随机信号模拟领域，具体地说，涉及一种非平稳随机过程高效模拟方法，所述的模拟方法提出了一种基于本征正交分解(POD)采用快速傅立叶变换(FFT)的谱表示模拟方法：首先对演化功率谱矩阵进行Cholesky分解，然后采用本征正交分解(POD)将分解出的时变谱表示为若干个时间函数与频率函数乘积之和，将时变功率谱进行解耦，最后利用FFT技术大幅度提高模拟的效率。采用本发明来模拟非平稳随机过程时，POD分解和重构的精度高，模拟出样本时程的统计相关函数与目标相关函数吻合度高；FFT使用使样本的模拟速度得到显著提高，本发明有效的解决了传统基于谱表示的非平稳随机过程模拟方法由于不能使用FFT而存在计算效率低下的问题，具备高效、精度高和易于使用的优点。

Description

一种非平稳随机过程高效模拟方法

技术领域

本发明属于随机信号模拟领域，具体地说，涉及一种非平稳随机过程高效模拟方法。

背景技术

在工程实践中，诸如极值风、地震动以及海浪等激励都呈现出明显的非平稳特征。而这些激励通常可被描述为零均值的非平稳随机过程，并且用Priestly所定义的演化功率谱(EPSD)来描述(Priestley M B.Evolutionary spectra andnon-stationary processes[J].Journal of Royal Statistical Society,1965,Series B,27:204-237.)。在采用时域方法进行非平稳随机响应分析时，基于演化功率谱的非平稳随机过程模拟是不可或缺的。

由于方法的精确性，基于谱表示法的非平稳随机过程模拟被广泛应用在工程领域，如地震动(钟菊芳，赵思烔，胡晓.基于时变功率谱的非平稳加速度时程合成研究[J].水利发电,2011,37(8):85-89.ZHONG Jufang,ZHAO Sitong,HUXiao.Method of nonstationary ground motion synthesizing based on evolutionarypower spectra[J].Water Power,2011,37(8):85-89.)和极值风(Huang G,Chen X,Liao H,Li M.Predicting tall building response to nonstationary winds using multiplewind speed samples[J].Wind and Structure,2013a,17(2):227-244.)的模拟等。然而，由于演化功率谱的时变性，该方法不能采用FFT技术来提高模拟效率(Deodatis G.Non-stationary stochastic vector processes:seismic ground motionapplications[J].Probabilistic Engineering Mechanics,1996b,11:149-168.)，从而限制了其在工程领域中的应用。为此，研究者通过对演化谱进行处理，提出了一些可以使用FFT的快速算法，如Li和Kareem基于随机分解概念提出的模拟方法和Huang提出的基于小波解耦的谱表示模拟方法(Huang G.An efficientsimulation approach for multivariate nonstationary process:hybrid of wavelet andspectral representation method[J].Probabilistic Engineering Mechanics,2014,37,74-83.)等，但这些方法存在精度相对不高、效率相对较低等缺点。

本征正交分解(POD)的理论基础是Karhunen-Loève(KL)展开，它被广泛应用于降维和数据的特征提取，如随机有限元分析(Ghanem R,Spanos P D.Stochasticfinite elements:A spectral approach[M].1991,Springer,Berlin.)和风压的压缩与重构(李元齐，沈祖炎.本征正交分解在曲面模型风场重构中的应用[J].同济大学学报(自然科学版),2006,34(1):22-26.LI Yuanqi,SONG Zuyan.Application of theproper orthogonal decomposition method to wind field reconstruction of modelswith curved surfaces[J].Journal of Tongji University(Natural Science),2006,34(1):22-26.)等。此外，不管是时域还是频域，在相关多点的平稳随机过程模拟中，POD的应用都是非常有效的(Chen X,Kareem A.Proper orthogonaldecomposition-based modeling,analysis,and simulation of dynamic wind loadeffects on structures[J].Journal of Wind Engineering and Industrial Aerodynamics,2005,131(4):1-18.)。而在非平稳随机过程模拟中，由于Fredholm积分很难处理，POD的应用则受到很大限制(Huang S P,Quek S T,Phoon K K.Convergence studyof the truncated Karhunen-Loeve expansion for simulation of stochastic processes[J].International Journal for Numerical Methods in Engineering,2001,52(9):1029-1043.)。

发明内容

鉴于此，本发明的目的是提供一种非平稳随机过程高效模拟方法，提出一种基于本征正交分解(POD)的采用快速傅立叶变换(FFT)的谱表示模拟方法，在该方法中首先对演化功率谱矩阵进行Cholesky分解，然后采用本征正交分解(POD)将分解出的时变谱表示为若干时间函数与频率函数乘积之和，即将时变功率谱进行解耦，最后利用FFT技术提高模拟的效率，采用本发明来模拟非平稳随机过程时，POD分解和重构的精度高，模拟出样本时程的统计相关函数与目标相关函数吻合度高，使FFT的使用使样本的模拟速度得到大幅度提高，有效的解决了传统基于谱表示的非平稳随机过程模拟方法由于不能使用FFT而存在计算效率低下的问题，具备高效、精度高和易于使用的优点。

为实现上述目的，本发明采用的技术方案为：一种非平稳随机过程高效模拟方法，其模拟方法的具体步骤如下：

1).获取非平稳随机过程演变功率谱密度EPSD的演化功率谱矩阵：

获取零均值n维向量过程x(t)＝[x₁(t),x₂(t),…,x_n(t)]^T的演化功率谱矩阵，如下式所示：

S(ω,t)＝[S_jk(ω,t)],j,k＝1,2,…,n                (1)

相干函数矩阵Γ(ω)满足下式：

Γ(ω)＝[γ_jk(ω)],j,k＝1,2,…,n              (2)

其中γ_jk(ω)是x_j(t)与x_k(t)之间的相干函数，x_j(t)与x_k(t)之间的互相关函数如下：

$R_{\mathrm{jk}}(t,t+\mathrm{\τ})={\∫}_{-\∞}^{\∞}\sqrt{S_{\mathrm{jj}}(\mathrm{\ω},t)S_{\mathrm{kk}}(\mathrm{\ω},t+\mathrm{\τ})}{\mathrm{\γ}}_{\mathrm{jk}}\left(\mathrm{\ω}\right)e^{\mathrm{i\ω\τ}}\mathrm{d\ω}---\left(3\right)$

2).对演化功率谱矩阵进行Cholesky分解，得到分解时变谱：

通过Cholesky分解将演化功率谱矩阵分解为上三角矩阵和下三角矩阵的乘积，得到下三角阵H(ω,t)为时变谱，如下式所示：

S(ω,t)＝H(ω,t)H^T*(ω,t)               (4)

其中*表示取复数共轭，H(ω,t)的对角元素为非负实数，其矩阵元素写成如下的复数形式，其中Im和Re分别表示虚部和实部；

$H_{\mathrm{jk}}(\mathrm{\ω},t)=\left|H_{\mathrm{jk}}(\mathrm{\ω},t)\right|e^{i{\mathrm{\θ}}_{\mathrm{jk}}\left(\mathrm{\ω}\right)}---\left(6\right)$

${\mathrm{\θ}}_{\mathrm{jk}}\left(\mathrm{\ω}\right)=\mathrm{ta}{n}^{-1}\left\{\frac{\mathrm{Im}\left[H_{\mathrm{jk}}(\mathrm{\ω},t)\right]}{\mathrm{Re}\left[H_{\mathrm{jk}}(\mathrm{\ω},t)\right]}\right\}---\left(7\right)$

3).将分解后的时变谱采用本征正交分解法(POD)分解：

将分解出的时变谱H_jk(ω,t)进一步表达为N_q(N_q≤8)个时间函数与频率函数乘积之和，如下所示：

$H_{\mathrm{jk}}(\mathrm{\ω},t)\≈\underset{q=1}{\overset{N_q}{\mathrm{\Σ}}}{a}_q\left(t\right){\mathrm{\Φ}}_q^{\mathrm{jk}}\left(\mathrm{\ω}\right),j,k=\mathrm{1,2},...,n;j\≥k---\left(8\right)$

其中a_q(t)是实的时间函数；为复数频率函数，如下式所示：

针对单变量随机过程和多变量随机过程求出a_q(t)和的值；

4).随机过程模拟：将非平稳随机过程模拟转化为平稳随机过程模拟：

将与H_jj(ω,t),j＝1,2,…,n对应的非平稳随机过程模拟转变为若干与有关的平稳随机过程模拟，假设满足下式：

又由于在方程(8)中，所有的频率函数是完全相干的，因此x_j(t)可用下式模拟：

$x_j\left(t\right)\≈\underset{m=1}{\overset{j}{\mathrm{\Σ}}}\underset{q=1}{\overset{N_q}{\mathrm{\Σ}}}{a}_q\left(t\right)y_q^{\mathrm{jm}}\left(t\right)---\left(12\right)$

从式(12)可以看出，可以看成是功率谱为的平稳随机过程模拟；

5).使用FFT高效模拟：

对转化后的平稳随机过程模拟，使用FFT技术来提高模拟效率。

作为优选，在步骤3)中将分解后的时变谱采用本征正交分解法(POD)分解时，对于单变量随机过程，对实数H₁₁(ω,t)作离散处理，采用通用的POD技术得出a_q(t)和的值。

作为优选，在步骤3)中将分解后的时变谱采用本征正交分解法(POD)分解时，对于多变量随机过程，采用以下处理步骤：

因为向量随机过程中的分量过程x_j(t)来自于同一类随机过程，因此所有S_jj(ω,t),j＝1,2,…,n具有相似的演化谱，从H_jj(ω,t),j＝1,2,…,n中分解出来的所有q阶时间函数都具有相似的形式，且所有的时间函数可用以下公式计算：

$a_q\left(t\right)\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{a}_q^{\mathrm{jj}}\left(t\right)---\left(13\right)$

则可以通过以下方式获得：在方程(8)两边同时乘于可得：

其中，

$Q_{q,r}={\∫}_0^{\∞}{a}_q\left(t\right)a_r\left(t\right)d_t,q,r=\mathrm{1,2},...,N_q---\left(15\right)$

$P_q^{\mathrm{jk}}\left(\mathrm{\ω}\right)={\∫}_0^{\∞}{H}_{\mathrm{jk}}(\mathrm{\ω},t)a_q\left(t\right)\mathrm{dt}---\left(16\right)$

当j＝k时，分别对H_jk(ω,t)的实部和虚部求解式(15)即可得出的实部和虚部，此时即求得a_q(t)和

本发明的有益效果为：

采用本发明来模拟非平稳随机过程时，POD分解和重构的精度高，模拟出样本时程的统计相关函数与目标相关函数吻合度高，使FFT的使用使样本的模拟速度得到大幅度提高，有效的解决了传统基于谱表示的非平稳随机过程模拟方法由于不能使用FFT而存在计算效率低下的问题，具备高效、精度高和易于使用的优点。

附图说明

图1为实例2的主坐标值和特征向量；

图2为实例2的近似谱与目标谱的对比；

图3为实例3模拟出的多点样本时程；

图4为实例3估计的相关函数对比(500个样本)。

具体实施方式

为了使本发明的目的、技术方案和有益效果更加清楚，下面将结合附图和实施例，对本发明的实施例进行详细的说明，以方便技术人员理解。

实施例1：具体模拟方法

一种非平稳随机过程高效模拟方法，其模拟方法的具体步骤如下：

1).获取非平稳随机过程演变功率谱密度EPSD的演化功率谱矩阵：

获取零均值n维向量过程x(t)＝[x₁(t),x₂(t),…,x_n(t)]^T的演化功率谱矩阵，如下式所示：

S(ω,t)＝[S_jk(ω,t)],j,k＝1,2,…,n                  (3)

相干函数矩阵Γ(ω)满足下式：

Γ(ω)＝[γ_jk(ω)],j,k＝1,2,…,n                 (4)

其中γ_jk(ω)是x_j(t)与x_k(t)之间的相干函数，x_j(t)与x_k(t)之间的互相关函数如下：

$R_{\mathrm{jk}}(t,t+\mathrm{\τ})={\∫}_{-\∞}^{\∞}\sqrt{S_{\mathrm{jj}}(\mathrm{\ω},t)S_{\mathrm{kk}}(\mathrm{\ω},t+\mathrm{\τ})}{\mathrm{\γ}}_{\mathrm{jk}}\left(\mathrm{\ω}\right)e^{\mathrm{i\ω\τ}}\mathrm{d\ω}---\left(3\right)$

2).对演化功率谱矩阵进行Cholesky分解，得到分解时变谱：

由于S(ω,t)在每一时刻t上，都是非负定的Hermitian矩阵，因此可通过Cholesky分解将演化功率谱矩阵分解为上三角矩阵和下三角矩阵的乘积，如下式所示：

S(ω,t)＝H(ω,t)H^T*(ω,t)                    (4)

其中H(ω,t)是下三角阵；*表示取复数共轭，H(ω,t)的对角元素为非负实数，非对角元素通常是复数，因此，其矩阵元素可写成如下的复数形式，其中Im和Re分别表示虚部和实部。

$H_{\mathrm{jk}}(\mathrm{\ω},t)=\left|H_{\mathrm{jk}}(\mathrm{\ω},t)\right|e^{i{\mathrm{\θ}}_{\mathrm{jk}}\left(\mathrm{\ω}\right)}---\left(6\right)$

${\mathrm{\θ}}_{\mathrm{jk}}\left(\mathrm{\ω}\right)=\mathrm{ta}{n}^{-1}\left\{\frac{\mathrm{Im}\left[H_{\mathrm{jk}}(\mathrm{\ω},t)\right]}{\mathrm{Re}\left[H_{\mathrm{jk}}(\mathrm{\ω},t)\right]}\right\}---\left(7\right)$

则向量过程中的每个元素过程x_j(t)可用如下公式模拟：

$x_j\left(t\right)=2\underset{m=1}{\overset{j}{\mathrm{\Σ}}}\underset{l=1}{\overset{N}{\mathrm{\Σ}}}\left|H_{\mathrm{jm}}({\mathrm{\ω}}_l,t)\right|\sqrt{\mathrm{\Δ\ω}}\cos [{\mathrm{\ω}}_{l}t-{\mathrm{\θ}}_{\mathrm{jm}}\left({\mathrm{\ω}}_l\right)+{\mathrm{\φ}}_{\mathrm{ml}}]---\left(17\right)$

式中Δω＝ω_u/N，ω_u是截止频率；ω_l＝lΔω；φ_1l,φ_2l,…,φ_nl为相互独立，且在[0,2π]服从均匀分布的随机相位角，由于H_jm(ω,t)是耦合的，因此式(17)不能使用FFT技术提高模拟效率。

3).将分解后的时变谱采用本征正交分解法(POD)分解：

将分解出的时变谱H_jk(ω,t)进一步表达为N_q(N_q≤8)个时间函数与频率函数乘积之和，如下所示：

$H_{\mathrm{jk}}(\mathrm{\ω},t)\≈\underset{q=1}{\overset{N_q}{\mathrm{\Σ}}}{a}_q\left(t\right){\mathrm{\Φ}}_q^{\mathrm{jk}}\left(\mathrm{\ω}\right),j,k=\mathrm{1,2},...,n;j\≥k---\left(8\right)$

其中a_q(t)是实的时间函数；为复数频率函数，如下式所示：

对于单变量随机过程，对实数H₁₁(ω,t)作离散处理，采用通用的POD技术得出a_q(t)和的值。

对于多变量随机过程，采用以下处理步骤：

因为向量随机过程中的分量过程x_j(t)来自于同一类随机过程，因此所有S_jj(ω,t),j＝1,2,…,n具有相似的演化谱，从H_jj(ω,t),j＝1,2,…,n中分解出来的所有q阶时间函数都具有相似的形式，且所有的时间函数可用以下公式计算：

$a_q\left(t\right)\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{a}_q^{\mathrm{jj}}\left(t\right)---\left(13\right)$

则可以通过以下方式获得：在方程(8)两边同时乘于可得：

其中，

$Q_{q,r}={\∫}_0^{\∞}{a}_q\left(t\right)a_r\left(t\right)d_t,q,r=\mathrm{1,2},...,N_q---\left(15\right)$

$P_q^{\mathrm{jk}}\left(\mathrm{\ω}\right)={\∫}_0^{\∞}{H}_{\mathrm{jk}}(\mathrm{\ω},t)a_q\left(t\right)\mathrm{dt}---\left(16\right)$

当j＝k时，分别对H_jk(ω,t)的实部和虚部求解式(15)即可得出的实部和虚部，此时即求得a_q(t)和

4).随机过程模拟：

将非平稳随机过程模拟转化为平稳随机过程模拟：

将与H_jj(ω,t),j＝1,2,…,n对应的非平稳随机过程模拟转变为若干与有关的平稳随机过程模拟，假设满足下式：

又由于在方程(8)中，所有的频率函数是完全相干的，因此x_j(t)可用下式模拟：

$x_j\left(t\right)\≈\underset{m=1}{\overset{j}{\mathrm{\Σ}}}\underset{q=1}{\overset{N_q}{\mathrm{\Σ}}}{a}_q\left(t\right)y_q^{\mathrm{jm}}\left(t\right)---\left(12\right)$

从式(12)可以看出，可以看成是功率谱为的平稳随机过程模拟；

5).使用FFT高效模拟：

对转化后的平稳随机过程模拟，使用FFT技术来提高模拟效率。

由上述模拟过程可知，如果忽略时变谱H_jk(ω,t)分解所需要的少量计算时间，本文提出的模拟方法只需要进行n(n+1)N_q/2次的FFT运算。而传统的谱表示模拟方法则需要进行n(n+1)/2次余弦项求和运算。

由于一次FFT计算的时间复杂度为O(Mlog₂M)(M为模拟时程的样本总数)，而直接求和的时间复杂度为O(M²)。因此，与常规的谱表示模拟方法相比，本文所提方法的计算效率得到了数十倍甚至百倍的提升。

实施例2：精度演示

基于POD的时变谱分解与重构，不失一般性，本文选取不可分离的演化功率谱进行算例分析，以验证采用POD进行时变谱分解的精确性，选取有代表性的单边功率谱密度函数如下式所示(Spanos P D,Kougioumtzoglou I A.Harmonicwavelets based statistical linearization for response evolutionary power spectrumdetermination[J].Probabilistic Engineering Mechanics,2012,27(1):57-68.)：

$S(f,t)=4\mathrm{\π}{(f/C)}^2{e}^{-\mathrm{Dt}}{t}^2{e}^{-{(f/C)}^{2_t}}---\left(18\right)$

其中，C＝2.5；D＝0.15。

对上述演化谱采用前述的POD分解方法，即可得出各阶模态的主坐标(时间函数)和特征向量(频率函数)。图1给出了前4阶模态的主坐标和特征向量。从图中可以看出，模态阶数越低，其时间函数包含的能量越大，且随着模态阶数的增大，能量衰减的很快。此外，在本次算例中，前6阶模态的特征值分别为12.207、3.787、0.582、0.140、0.018和0.003。因此，采用前几阶模态就可包含绝大部分的能量。由此看见，采用POD进行时变谱的分解和重构具有很明显的物理意义。更重要的是，POD是基于数据的最优滤波器，因此该方法在时变谱的分解中非常有效。图2为采用前6阶模态重构出的近似功率谱与目标谱的对比图，从图中可以看出，两者之间几乎没有区别。

实施例3：效果分析

多点非平稳随机过程的模拟实例

为了更好地说明本文所提方法的优越性，本文选取了(Li Y,Kareem A.Simulation of multivariate nonstationary random processes by FFT[J].Journal ofEngineering Mechanics,1991,117:1037-1058.)中使用的算例进行多点非平稳随机过程模拟，并将本文计算结果与该文献计算结果进行对比。在该算例中，由于土壤的液化作用，地震动的演化功率谱密度在7s以后存在突变，即由宽带的频谱变成集中在低频的窄带频谱。

设{x(t)}＝[x₁(t),x₂(t)]^T，且x₁(t)和x₂(t)的演化功率谱密度都服从下式：

$S(f,t)=A_1^2\left(t\right)K(f,t)---\left(19\right)$

其中，

$A_1\left(t\right)=\frac{\exp (-\mathrm{ct})-\exp (-\mathrm{dt})}{\max [\exp (-\mathrm{ct})-\exp (-\mathrm{dt})]}---\left(20\right)$

且c＝0.25，d＝0.5；K(f,t)满足下式：

$K(f,t)=4\mathrm{\π}{S}_0\frac{f_g^4\left(t\right)+4{\mathrm{\ζ}}_g^2\left(t\right)f_g^2\left(t\right)f^2}{{(f^2-f_g^2\left(t\right))}^2+4{\mathrm{\ζ}}_g^2\left(t\right)f_g^2\left(t\right)f^2}---\left(21\right)$

其中， $f_g\left(t\right)=\left\{\begin{array}{cc}2.476\mathrm{Hz}& \mathrm{for}0\&le;t\&le;4.5s\\ 4.1316t^{\&prime;3}& -6.474t^{\&prime;2}+2.467\mathrm{Hz}\\ & \mathrm{for}4.5s<t<5.5s\\ 0.3183\mathrm{Hz}& \mathrm{for}5.5s\&le;t\end{array}\right.{\mathrm{\&zeta;}}_g\left(t\right)\left\{\begin{array}{cc}0.64& \mathrm{for}0\&le;t\&le;4.5s\\ 1.250t^{\&prime;3}& -1.875t^{\&prime;2}+0.64\\ & \mathrm{for}4.5s<t<5.5s\\ 0.015& \mathrm{for}5.5s\&le;t\end{array}\right.---\left(22\right)$ 在式(22)中，t'＝t-4.5。x₁(t)和x₂(t)的相关性如下式所示^[14]：

$\mathrm{Co}{h}_{12}^0\left(f\right)=\mathrm{\&rho;}({\mathrm{\&upsi;}}_{12},f)\exp (-i2\mathrm{\&pi;f}\frac{{\mathrm{\&upsi;}}_{12}}{V_{12}})---\left(23\right)$

$\mathrm{\&rho;}({\mathrm{\&upsi;}}_{12},f)=\mathrm{\&beta;exp}[-\frac{2{\mathrm{\&upsi;}}_{12}}{\mathrm{\&alpha;\&theta;}\left(f\right)}(1-\mathrm{\&beta;}-\mathrm{\&alpha;\&beta;})]+(1-\mathrm{\&beta;})\exp [-\frac{2{\mathrm{\&upsi;}}_{12}}{\mathrm{\&theta;}\left(f\right)}(1-\mathrm{\&beta;}-\mathrm{\&alpha;\&beta;})]---\left(24\right)$

θ(f)＝θ₀[1+(f/f₀)^b]^-1/2                     (25)

其中υ₁₂表示两点间的距离，且υ₁₂＝150m；V₁₂为地震波沿两点间直线方向的传播速度，且V₁₂＝2000m/s；α＝0.147；β＝0.736；θ₀＝5210m/s；f₀＝1.09Hz。

在随机过程模拟中，选取的参数如下：截止频率f_u＝10Hz；频域上的样本个数N＝1024；时间间隔Δt＝0.05s；样本个数M＝2048。图3为模拟出的非平稳样本时程，图4则为时间延迟在τ＝0s，0.1s以及0.2s上的估计相关函数与目标相关函数的对比。由图中可以看出，两者之间吻合很好。值得说明的是，与Li和Kareem^[5]的计算结果相比，本文所提方法具有很大的优势。这是因为，本文方法可以捕捉出目标相关函数在5-6s之间的峰值，且由样本估计出的统计相关函数更加平滑。此外，本文POD分解使用的模态阶数N_q＝8，而Li和Kareem则需要用15阶多项式去拟合目标值。

综上所述，本发明提供一种非平稳随机过程高效模拟方法，提出一种基于本征正交分解(POD)的采用快速傅立叶变换(FFT)的谱表示模拟方法，在该方法中首先对演化功率谱矩阵进行Cholesky分解，然后采用本征正交分解(POD)将分解出的时变谱表示为若干时间函数与频率函数乘积之和，即将时变功率谱进行解耦，最后利用FFT技术提高模拟的效率，采用本发明来模拟非平稳随机过程时，POD分解和重构的精度高，模拟出样本时程的统计相关函数与目标相关函数吻合度高，使FFT的使用使样本的模拟速度得到大幅度提高，有效的解决了传统基于谱表示的非平稳随机过程模拟方法由于不能使用FFT而存在计算效率低下的问题，具备高效、精度高和易于使用的优点。

需要说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (3)

1.一种非平稳随机过程高效模拟方法，其特征在于：所述的模拟方法的具体步骤如下：

1).获取非平稳随机过程演变功率谱密度EPSD的演化功率谱矩阵：

获取零均值n维向量过程x(t)＝[x₁(t),x₂(t),…,x_n(t)]^T的演化功率谱矩阵，如下式所示：

S(ω,t)＝[S_jk(ω,t)],j,k＝1,2,…,n  (1)

相干函数矩阵Γ(ω)满足下式：

Γ(ω)＝[γ_jk(ω)],j,k＝1,2,…,n  (2)

其中γ_jk(ω)是x_j(t)与x_k(t)之间的相干函数，x_j(t)与x_k(t)之间的互相关函数如下：

$R_{\mathrm{jk}}(t,t+\mathrm{\&tau;})={\&Integral;}_{-\&infin;}^{\&infin;}\sqrt{S_{\mathrm{jj}}(\mathrm{\&omega;},t+\mathrm{\&tau;})}{\mathrm{\&gamma;}}_{\mathrm{jk}}\left(\mathrm{\&omega;}\right)e^{\mathrm{i\&omega;\&tau;}}\mathrm{d\&omega;}---\left(3\right)$

2).对演化功率谱矩阵进行Cholesky分解，得到分解时变谱：

通过Cholesky分解将演化功率谱矩阵分解为上三角矩阵和下三角矩阵的乘积，得到下三角阵H(ω,t)为时变谱，如下式所示：

S(ω,t)＝H(ω,t)H^T*(ω,t)  (4)

其中*表示取复数共轭，H(ω,t)的对角元素为非负实数，其矩阵元素写成如下的复数形式，其中Im和Re分别表示虚部和实部；

$H_{\mathrm{jk}}(\mathrm{\&omega;},t)=\left|H_{\mathrm{jk}}(\mathrm{\&omega;},t)\right|e^{{\mathrm{j\&theta;}}_{\mathrm{jk}}\left(\mathrm{\&omega;}\right)}---\left(6\right)$

${\mathrm{\&theta;}}_{\mathrm{jk}}\left(\mathrm{\&omega;}\right)={\tan }^{-1}\left\{\frac{\mathrm{Im}\left[H_{\mathrm{jk}}(\mathrm{\&omega;},t)\right]}{\mathrm{Re}\left[H_{\mathrm{jk}}(\mathrm{\&omega;},t)\right]}\right\}---\left(7\right)$

3).将分解后的时变谱采用本征正交分解法(POD)分解：

将分解出的时变谱H_jk(ω,t)进一步表达为N_q(N_q≤8)个时间函数与频率函数乘积之和，如下所示：

$H_{\mathrm{jk}}(\mathrm{\&omega;},t)\&ap;\underset{q=1}{\overset{N_q}{\mathrm{\&Sigma;}}}{a}_q\left(t\right){\mathrm{\&Phi;}}_q^{\mathrm{jk}}\left(\mathrm{\&omega;}\right),j,k=\mathrm{1,2},...,n;j\&le;k---\left(8\right)$

其中a_q(t)是实的时间函数；为复数频率函数，如下式所示：

${\mathrm{\&Phi;}}_q^{\mathrm{jk}}\left(\mathrm{\&omega;}\right)={\mathrm{\&Phi;}}_q^{\mathrm{jk}}\left(\mathrm{\&omega;}\right)e^{{\mathrm{i\&theta;}}_q^{\mathrm{jk}}}\left(\mathrm{\&omega;}\right)---\left(9\right)$

$l_q^{\mathrm{jk}}\left(\mathrm{\&omega;}\right)={\tan }^{-1}\left\{\frac{\mathrm{Im}\left[{\mathrm{\&Phi;}}_q^{\mathrm{jk}}\left(\mathrm{\&omega;}\right)\right]}{\mathrm{Re}\left[{\mathrm{\&Phi;}}_q^{\mathrm{jk}}\left(\mathrm{\&omega;}\right)\right]}\right\}---\left(10\right)$

针对单变量随机过程和多变量随机过程求出a_q(t)和的值；

4).随机过程模拟：将非平稳随机过程模拟转化为平稳随机过程模拟：

将与H_jj(ω,t),j＝1,2,…,n对应的非平稳随机过程模拟转变为若干与有关的平稳随机过程模拟，假设满足下式：

又由于在方程(8)中，所有的频率函数是完全相干的，因此x_j(t)可用下式模拟：

$x_j\left(t\right)\&ap;\underset{m=1}{\overset{j}{\mathrm{\&Sigma;}}}\underset{q=1}{\overset{N_q}{\mathrm{\&Sigma;}}}{a}_q\left(t\right)y_q^{\mathrm{jm}}\left(t\right)---\left(12\right)$

从式(12)

可以看出，可以看成是功率谱为的平稳随机过程模拟；

5).使用FFT高效模拟：

对转化后的平稳随机过程模拟，使用FFT技术来提高模拟效率。

2.根据权利要求1所述的一种非平稳随机过程高效模拟方法，其特征在于：在步骤3)中将分解后的时变谱采用本征正交分解法(POD)分解时，对于单变量随机过程，对实数H₁₁(ω,t)作离散处理，采用通用的POD技术得出a_q(t)和的值。

3.根据权利要求1所述的一种非平稳随机过程高效模拟方法，其特征在于：在步骤3)中将分解后的时变谱采用本征正交分解法(POD)分解时，对于多变量随机过程，采用以下处理步骤：

因为向量随机过程中的分量过程x_j(t)来自于同一类随机过程，因此所有S_jj(ω,t),j＝1,2,…,n具有相似的演化谱，从H_jj(ω,t),j＝1,2,…,n中分解出来的所有q阶时间函数都具有相似的形式，且所有的时间函数可用以下公式计算：

$a_q\left(t\right)=\frac{1}{n}\underset{j=1}{\overset{n}{\mathrm{\&Sigma;}}}{a}_q^{\mathrm{jj}}\left(t\right)---\left(13\right)$

则可以通过以下方式获得：在方程(8)两边同时乘于可得：

其中，

$Q_{q,r}={\&Integral;}_0^{\&infin;}{a}_q\left(t\right)a_r\left(t\right)\mathrm{dt},q,r=\mathrm{1,2},...,N_q---\left(15\right)$

$P_q^{\mathrm{jk}}\left(\mathrm{\&omega;}\right)={\&Integral;}_0^{\&infin;}{H}_{\mathrm{jk}}(\mathrm{\&omega;},t)a_q\left(t\right)\mathrm{dt}---\left(16\right)$

当j＝k时，分别对H_jk(ω,t)的实部和虚部求解式(15)即可得出的实部和虚部，此时即求得a_q(t)和