CN107247687A - 一种基于特征正交分解的非平稳随机过程快速模拟方法 - Google Patents

一种基于特征正交分解的非平稳随机过程快速模拟方法 Download PDF

Info

Publication number
CN107247687A
CN107247687A CN201710351315.4A CN201710351315A CN107247687A CN 107247687 A CN107247687 A CN 107247687A CN 201710351315 A CN201710351315 A CN 201710351315A CN 107247687 A CN107247687 A CN 107247687A
Authority
CN
China
Prior art keywords
mrow
msub
msubsup
munderover
omega
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201710351315.4A
Other languages
English (en)
Other versions
CN107247687B (zh
Inventor
赵宁
黄国庆
陆成文
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest Jiaotong University
Original Assignee
Southwest Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Southwest Jiaotong University filed Critical Southwest Jiaotong University
Priority to CN201710351315.4A priority Critical patent/CN107247687B/zh
Publication of CN107247687A publication Critical patent/CN107247687A/zh
Application granted granted Critical
Publication of CN107247687B publication Critical patent/CN107247687B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G06F17/156—Correlation function computation including computation of convolution operations using a domain transform, e.g. Fourier transform, polynomial transform, number theoretic transform

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Computational Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Discrete Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)

Abstract

本发明公开了一种基于特征正交分解的非平稳随机过程快速模拟方法，对于空间任意分布的、时变相干的多点非平稳过程，首先对EPSD矩阵进行三角分解，然后用POD将三角矩阵的各个时频耦合的元素分别表示为若干个时间函数和频率函数乘积之和，最后采用FFT技术进行快速模拟；对于空间任意分布的、时不变相干的多点非平稳过程，先对相干矩阵进行三角分解，再使用POD将自谱组成的对角矩阵的时频耦合的元素，分别表示为若干个时间函数和频率函数乘积之和，最后优化模拟公式，用极少次的FFT技术进行高效模拟。本发明采用POD拟合精度高，模拟出的样本准确；三角分解及FFT技术的优化使用使模拟效率大大提升；解决了目前非平稳模拟方法存在的应用局限及效率低下等问题。

Description

一种基于特征正交分解的非平稳随机过程快速模拟方法
技术领域
本发明涉及随机信号模拟技术领域，具体为一种基于特征正交分解的非平稳随机过程快速模拟方法。
背景技术
地震和瞬时极端风等毁灭性的激励表现出非平稳特性，通常用非平稳过程来描述。当涉及到非线性、系统随机性、参数激励和其他某些随机问题时，非平稳过程的蒙特卡洛模拟(MCS)是结构响应分析中的一项重要任务。此外，MCS是评价其他随机方法准确性的基准。
基于演化功率谱密度函数(EPSD Evolutionary Power Spectral Density)的谱表示方法，由于其准确性和简便性被广泛应用于随机过程的模拟中。然而，EPSD的时变特性使模拟出现了两个困难。一方面，功率谱矩阵在每个时刻和频率都要进行一次三角分解。另一方面，由于不能使用快速傅里叶变换(FFT)，非平稳随机过程的模拟效率较为低下。
(Li,Y.,and Kareem,A.Simulation of multivariate nonstationary randomprocesses by FFT[J].Journal of Engineering Mechanics,1991,117,1037-58.)使用了三角函数和勒让德多项式函数来拟合分解后的EPSD矩阵，建立了使用FFT的谱表示方法。(Huang,G.An efficient simulation approach for multivariate nonstationaryprocess:Hybrid of wavelet and spectral representation method[J].Probabilistic.Engineering Mechanics.,2014,37,74-83.)使用小波解耦时变EPSD矩阵，发展了一种基于快速傅里叶变换的谱表示方法。(Huang,G.Application of ProperOrthogonal Decomposition in Fast Fourier Transform—Assisted MultivariateNonstationary Process Simulation[J].Journal of Engineering Mechanics 2015,141.7:04015015.)采用特征正交分解(POD Proper Orthogonal Decomposition)来分离分解后的EPSD矩阵，并利用快速傅里叶变换来提高模拟速度。这里的基于POD的分解是基于数据的，最优的并且是易于工程应用的。(Peng,L.,Huang,G.,Chen,X.,and Kareem,A.Simulation of multivariate nonstationary random process:A hybrid stochasticwave and proper orthogonal decomposition approach[J].Journal of EngineeringMechanics(2017accepted))提出了一种新的基于随机波和POD的混合方法，用于模拟沿直线分布的多点非平稳过程。其中二维FFT的使用大大地提高了模拟的效率。总的来说，在非平稳模拟方面，虽然已经取得了很大的进步，但是仍然有必要提高模拟的效率和适用性，特别是对于有许多模拟点的现代庞大的结构。
此外，一些非平稳激励的相干函数具有明显的时变特征，如下击暴流风。因此，有必要把常规的模拟方法拓展到时变相干的非平稳过程。
发明内容
针对上述问题，本发明的目的在于提供一种可有效地解决目前各种非平稳模拟方法存在的应用局限以及效率低下等问题的基于特征正交分解的非平稳随机过程快速模拟方法，采用该方法模拟出的样本准确，模拟效率高的。技术方案如下：
一种任意非平稳随机过程快速模拟方法，包括以下步骤：
步骤1：获取时变相干的多点非平稳随机过程的目标功率谱和相干函数：
获取零均值n维向量过程x(t)＝[x1(t),x2(t),…,xn(t)]T的演化功率谱矩阵，如下式所示：
其中，ω是圆频率；t是时间；Sjj(ω,t)是随机过程xj(t)的自谱，Sjk(ω,t)是两个随机过程xj(t)和xk(t)之间的互谱，通过下式确定：
其中：γjk(ω,t)是xj(t)和xk(t)的时变的相干函数，由相干函数组成的相干矩阵Γ(ω,t)由下式给出：
步骤2：对演化功率谱矩阵进行三角分解，得到分解时变谱，进而得到谱表示模拟公式：
通过Cholesky分解，演化功率谱矩阵分解为：
其中，H(ω,t)为下三角矩阵形，表示为：
任意两个随机过程之间的演化功率谱Sjk(ω,t)表达为：
它们的相关函数表达为：
那么谱表示模拟公式为：
其中，△ω＝ωu/N是频率分辨率，ωu是截止频率，N是离散频率数，ωl＝l△ω，φkl是在区间[0,2π]均匀分布的独立随机相位角；
步骤3：将得到的三角矩阵各个时变谱作特征正交分解：
将三角分解得到的时变谱Hjk(ω,t)进行特征正交分解后，再近似表达为多个时间函数和频率函数乘积之和：
其中：是频率相关矩阵Rjk的第q个特征向量；是通过计算的第q个主坐标；是包含绝大部分能量的有效项数；
步骤4：使用FFT技术模拟非平稳过程：
利用特征正交分解拟合式，将所述谱表示模拟公式，即式(8)转换为FFT直接可用的形式，如下所示：
再采用FFT技术快速模拟。
进一步的，所述步骤3中，对时变谱Hjk(ω,t)采用特征正交分解的具体过程下：
若使特征正交分解可使用，则Hjk(ω,t)需在时频域上离散化，该连续的二元函数离散为一个N维列向量：
其中，表示第l个频率点处的时间序列，通过特征正交分解，找到一组最优正交基为在这组正交基上，Hjk(p△t)的投影最大；这组基向量通过以下特征向量分解来获得：
其中，
其为时间平均的频率相关矩阵；是第q个特征向量，为频率的函数，是第q个特征值；当特征值按降序排列时，低阶特征值包含更多能量，则Hjk(p△t)近似表示为：
其中，是第q个主坐标或时间函数，用确定。
更进一步的，所述步骤4中，采用FFT技术对谱表示模拟公式进行快速计算，其具体处理如下：
为了便于FFT应用于模拟，模拟公式改写为如下形式：
其中，M是时间离散点数；为：
时间和频率分辨率满足：
M△t＝2π/△ω (17)
为了避免混叠，时间步长也需要满足：
△t≤2π/(2ωu) (18)
也即
M≥2N (19)
进一步改写为：
为充分利用FFT技术，时间的离散数满足M＝2u，u为正整数。
一种时不变相干的非平稳随机过程快速模拟方法，包括以下步骤：
步骤A：获取时不变相干的非平稳随机过程的目标功率谱和相干函数：
时不变相干的非平稳过程的演化功率谱矩阵表示为：
S(ω,t)＝D(ω,t)Γ(ω)DT(ω,t) (21)
其中
其中，ω是圆频率，γjk(ω)是随机过程xj(t)和xk(t)之间的时不变的相干函数，它们之间的相关函数表示为：
其中，Sjj(ω,t)为随机过程xj(t)的自谱，Skk(ω,t+t)为随机过程xk(t)的自谱；
步骤B：对相干矩阵进行三角分解，得到谱表示模拟公式：
相干矩阵Γ(ω)为正定的Hermitian矩阵，分解为：
其中，B(ω)为下三角矩阵，表示为：
那么谱表示模拟公式为：
其中，△ω＝ωu/N是频率分辨率，ωu是截止频率，N是离散频率数，ωl＝l△ω，φkl是在区间[0,2π]均匀分布的独立随机相位角；
步骤C：对关于自谱的元素作特征正交分解，得到谱表示模拟公式：
将元素近似表达为多个时间函数和频率函数乘积之和：
其中：是频率相关矩阵Rjj的第q个特征向量；是通过计算的第q个主坐标；是包含绝大部分能量的有效项数；
步骤D：优化模拟公式，并采用FFT技术模拟非平稳过程：
利用特征正交分解拟合式，将所述谱表示模拟公式，即式(27)转换为FFT直接可用的形式，然后交换公式求和的顺序，减少FFT的执行次数，如下所示：
其中，△ω＝ωu/N是频率分辨率，ωu是截止频率，N是离散频率数，ωl＝l△ω；
再采用FFT技术快速模拟非平稳过程时程样本。
本发明的有益效果是：本发明的方法能够模拟任意的非平稳过程，由于POD拟合的精度高，所以模拟出的样本准确；三角分解及FFT技术的优化使用使模拟效率大大提升；有效地解决了目前各种非平稳模拟方法存在的应用局限以及效率低下等问题，具有广阔的应用前景。
附图说明
图1a为实例1中的模拟点3和模拟点5分解后的时变谱的主坐标
图1b为实例1中的模拟点3和模拟点5分解后的特征向量。
图2a为实例1中模拟点3采用两种模拟方法拟合的时变谱与目标值的对比。
图2b为实例1中模拟点5采用两种模拟方法拟合的时变谱与目标值的对比。
图3为实例1中各模拟点的样本时程。
图4为实例1中模拟点3自相关函数与目标值的对比。
图5为实例1中模拟点5的自相关函数与目标值的对比。
图6为实例1中模拟点3模拟点5的互相关函数与目标值的对比。
具体实施方式
下面结合附图和具体实施方式对本发明作进一步详细的说明。本实施例示出一种基于特征正交分解(POD)采用快速傅立叶变换(FFT)的谱表示模拟方法：对于空间任意分布的、时变相干的多点非平稳过程，首先对演化功率谱矩阵进行三角分解，然后使用特征正交分解(POD)将三角矩阵的各个时频耦合的元素分别表示为若干个时间函数和频率函数乘积之和，最后采用FFT技术进行快速模拟。特别地，当该方法应用于空间任意分布的、时不变相干的多点非平稳过程时，首先对相干矩阵进行三角分解，然后使用POD将自谱组成的对角矩阵的时频耦合的元素，分别表示为若干个时间函数和频率函数乘积之和，最后优化模拟公式，用极少次的FFT技术进行高效模拟。
一种任意的非平稳随机过程快速模拟方法，具体步骤如下：
步骤1：获取多点非平稳随机过程的目标功率谱和相干函数：
获取零均值n维向量过程x(t)＝[x1(t),x2(t),…,xn(t)]T的演化功率谱(EPSD)矩阵，如下式所示：
其中：ω是圆频率；Sjj(ω,t)是随机过程xj(t)的自谱，Sjk(ω,t)是两个随机过程xj(t)和xk(t)之间的互谱，可通过下式确定：
其中：γjk(ω,t)是xj(t)和xk(t)的时变的相干函数。由相干函数组成的相干矩阵Γ(ω,t)由下式给出：
然后，对频率和时间离散化处理，进而获得在离散时频域上的演化功率谱矩阵。
步骤2：在离散的时频点上对演化功率谱矩阵进行三角分解，得到离散的分解时变谱。
通过Cholesky分解，EPSD矩阵可以作出以下分解：
其中：H(ω,t)是一个下三角矩阵形如：
它的对角元素一般是非负实数；而非对角元素通常是复数。
任意两个随机过程之间的演化功率谱Sjk(ω,t)表达为：
它们的相关函数表达为：
那么谱表示模拟公式为：
其中，△ω＝ωu/N是频率分辨率，ωu是截止频率，N是离散频率数，ωl＝l△ω，φkl是在区间[0,2π]均匀分布的独立随机相位角。
步骤3：将得到的三角矩阵各个离散的时变谱作特征正交分解：
将分解后的时变谱采用特征正交分解(POD)的具体处理如下：
为了使POD能够使用，Hjk(ω,t)应该在时频域上离散化。该连续的二元函数可以离散为一个N维列向量：
其中表示第l个频率点处的时间序列。通过POD，可以找到一组最优正交基为在这组正交基上，Hjk(p△t)的投影是最大的。这组基向量可以通过以下特征向量分解来获得，
其中
它可以被解释为时间平均的频率相关矩阵；是第q个特征向量，它也是一个频率的函数，是第q个特征值。当特征值按降序排列时，低阶特征值包含更多能量。那么，
Hjk(p△t)可以近似表示为：
其中，是第q个主坐标(或时间函数)，可以用确定。
将分解得到的时变谱Hjk(ω,t)进一步近似表达为多个时间序列和频率序列乘积之和，如下所示：
其中是频率相关矩阵Rjk的第q个特征向量；是通过计算的第q个主坐标；是选取的包含绝大部分能量的有效项数。
利用POD拟合式(42)，将谱表示模拟公式(37)转换为FFT直接可用的形式：
步骤4：使用FFT技术高效地模拟非平稳过程：
式(43)可进一步写为离散化的形式：
其中M是时间离散点数；△t是时间分辨率，时间分辨率△t和频率分辨率△ω满足：
M△t＝2π/△ω (45)
为了避免混叠，时间步长也需要满足：
△t≤2π/(2ωu) (46)
也即
M≥2N (47)
进一步改写为：
其中，满足M＝2u，u是一个正整数；为：
其中△ω＝ωu/N是频率分辨率，ωu是截止频率，N是离散频率数。然后使用FFT技术快速模拟出非平稳时程样本。进一步，通过下式确定目标相关函数：
然后，将模拟的样本计算的相关函数与目标相关函数进行对比，验证模拟的精度。
值得指出的是，经典的谱表示模拟公式将分解后的时变谱表示为极坐标的形式，如下：
其中|·|表示复数的模；θjk(ω,t)是Hjk(ω,t)的相位角，通过下式确定：
θjk(ω,t)＝tan-1{Im[Hjk(ω,t)]/Re[Hjk(ω,t)]} (52)
其中：Im和Re分别表示虚部和实部。但是这样求得的相位角只在(-π/2,π/2)的角度范围，而实际的相位角应在(-π,π)的范围内。但是，本发明所提方法避免了这样的问题。
对于时不变相干的非平稳过程的模拟，其具体实施步骤为：
步骤A：获取时不变非平稳随机过程的目标功率谱和相干函数：
时不变相干的非平稳过程的演化功率谱矩阵可以表示为：
S(ω,t)＝D(ω,t)Γ(ω)DT(ω,t) (53)
其中
其中：γjk(ω)是xj(t)和xk(t)之间的时不变的相干函数。然后，对频率和时间离散化处理，进而获得在离散时频域上的功率谱和相干函数。
它们之间的相关函数表示为：
其中，Sjj(ω,t)为随机过程xj(t)的自谱，Skk(ω,t+t)为随机过程xk(t)的自谱。
步骤B：对离散频域上的相干矩阵进行三角分解：
相干矩阵Γ(ω)也是一个正定的Hermitian矩阵，它可以被分解为
其中B(ω)是一个下三角矩阵，可以表示为
那么谱表示模拟公式为：
其中，△ω＝ωu/N是频率分辨率，ωu是截止频率，N是离散频率数，ωl＝l△ω，φkl是在区间[0,2π]均匀分布的独立随机相位角。
步骤C：对关于自谱的元素作特征正交分解：
将元素近似表达为多个时间序列和频率序列乘积之和，如下所示：
其中：是频率相关矩阵Rjj的第q个特征向量；是通过计算的第q个主坐标；是包含绝大部分能量的有效项数。
利用上述POD拟合式(60)，将谱表示模拟公式(59)转换为FFT直接可用的形式：
步骤D：优化模拟公式，并使用FFT技术高效地模拟非平稳过程：
模拟公式(61)可进一步写为离散化的形式，然后交换公式求和的顺序，减少FFT的执行次数，如下所示：
其中，
最后使用FFT技术快速模拟非平稳过程的时程样本。进一步目标相关函数可确定为：
然后，将模拟的样本计算的相关函数与目标相关函数进行对比，验证模拟的精度。
下面通过具体的实例对本发明以及本发明效果进行进一步的说明。
2015年,(Huang,Guoqing.Application of Proper Orthogonal Decompositionin Fast Fourier Transform—Assisted Multivariate Nonstationary ProcessSimulation[J].Journal of Engineering Mechanics 141.7(2015):04015015.)将POD应用于非平稳随机过程的谱表示方法，提出了一种可以采用FFT的非平稳随机过程快速模拟方法。以下通过两个实例分别与Huang的方法进行对比，验证本发明方法的精度和效率。
实例1：模拟精度验证
本例将讨论沿高层建筑下击暴流风场的模拟。首先，通过在0.9、2.4、4、10、116、158和200m高度处(编号1-7)实测的下击暴流风速数据，估计出该多点非平稳过程的演化功率谱和相干函数。从相干函数的特性发现它是时变的，然后按照以上时变相干的非平稳过程的模拟步骤进行样本生成。
图1a和图1b分别显示了点3和点5对应的分解后的时变谱H33(f,t)and H55(f,t)的前两阶主坐标和特征向量。可以发现：相应于较小特征值的主坐标(时间函数)包含了较大的能量。图2a和图2b比较了所提方法拟合的分解后的谱与相应的目标值。其中，使用POD时，仅仅采用前5阶就保留了99％的能量。因此，该方法在匹配分解后的谱方面具有优秀的表现。为了便于比较，由Huang的方法拟合的分解后的谱也绘制于图2a和图2b中，这种方法使用了前32阶进行POD拟合。从图中可以看到，这种方法在拟合分解后的谱时可能会带来较大误差。
在模拟中，截止频率被取为0.5Hz、频率增量为0.00049Hz。生成的时程样本的时间分辨率为1s，持续时间为2048s。图3显示了模拟的7个点的时程样本。然后用10000个模拟样本计算自相关/互相关函数。图4、图5、图6对比了在0、4、8秒三个时间差上，估算的相关函数与目标值。可以看出，本发明具有很高的模拟精度。
实例2：模拟效率对比
一座主跨为2000米，边跨为1000米的大跨度悬索桥沿桥面的风场被研究，为了简单起见，只考虑垂直于桥轴线方向的风速分量。边跨和主跨上两个相邻的模拟点之间的距离分别取为10米和20米。因此，共301个模拟点的风场将会被模拟。该风场特性根据实测的台风记录构造。同前面的例子一样，从实测的数据估计出台风的功率谱。然后，假定桥梁中跨模拟点的自谱等于该台风谱，而边跨各点的自谱取为0.8倍的台风谱。此外，相干函数采用了Davenport的指数函数模型：
其中：λ＝7是衰减因子，ξ是桥面上两点间距，U0＝40m/s是桥面高度处的平均风速。
在模拟中，截止频率取为0.5Hz，频率增量为0.000244Hz。生成的时程样本的时间分辨率为1s、持续时间为4096s。当要确保99％的能量时，POD拟合需要的近似项的数量为基于本算例我们对比了它和黄的方法的在这三方面的计算时间，如表1所示。其中，黄的方法中POD的计算时间包括了求解线性联立方程所用的时间。表中时间比是指黄的方法的计算时间与所提方法的时间的比值。为了不失一般性，对各点自谱都使用了POD。从表中可以发现，总的模拟效率提高了645倍。特别地，如果只对两个不同自谱应用POD，那么总的效率可提升3600倍。
表1实例2中两种方法模拟效率的对比(s)
计算过程 本发明方法 Huang方法 时间比
Cholesky分解 326 1173600 3600
POD 1592 67364 42
FFT 7 57 8
总时间 1925 1241021 645
此外，通过进一步计算，当有50、175、301个模拟点时，总的时间比分别是117、388、645。可以发现，模拟点数越多，效率提升越明显。

Claims (4)

1.一种任意非平稳随机过程快速模拟方法，其特征在于，包括以下步骤：
步骤1：获取时变相干的多点非平稳随机过程的目标功率谱和相干函数：
获取零均值n维向量过程x(t)＝[x1(t),x2(t),…,xn(t)]T的演化功率谱矩阵，如下式所示：
其中，ω是圆频率；t是时间；Sjj(ω,t)是随机过程xj(t)的自谱，Sjk(ω,t)是两个随机过程xj(t)和xk(t)之间的互谱，通过下式确定：
<mrow> <msub> <mi>S</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msqrt> <mrow> <msub> <mi>S</mi> <mrow> <mi>j</mi> <mi>j</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <msub> <mi>S</mi> <mrow> <mi>k</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </msqrt> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
其中：γjk(ω,t)是xj(t)和xk(t)的时变的相干函数，由相干函数组成的相干矩阵Γ(ω,t)由下式给出：
步骤2：对演化功率谱矩阵进行三角分解，得到分解时变谱，进而得到谱表示模拟公式：通过Cholesky分解，演化功率谱矩阵分解为：
S(ω,t)＝H(ω,t)HT*(ω,t) (4)
其中，H(ω,t)为下三角矩阵形，表示为：
任意两个随机过程之间的演化功率谱Sjk(ω,t)表达为：
<mrow> <msub> <mi>S</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>H</mi> <mrow> <mi>j</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>H</mi> <mrow> <mi>k</mi> <mi>m</mi> </mrow> <mo>*</mo> </msubsup> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>6</mn> <mo>)</mo> </mrow> </mrow>
它们的相关函数表达为：
<mrow> <msub> <mi>R</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>t</mi> <mo>+</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mo>&amp;Integral;</mo> <mrow> <mo>-</mo> <mi>&amp;infin;</mi> </mrow> <mi>&amp;infin;</mi> </msubsup> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>n</mi> </munderover> <msub> <mi>H</mi> <mrow> <mi>j</mi> <mi>m</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>H</mi> <mrow> <mi>k</mi> <mi>n</mi> </mrow> <mo>*</mo> </msubsup> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>+</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mi>i</mi> <mi>&amp;omega;</mi> <mi>&amp;tau;</mi> </mrow> </msup> <mi>d</mi> <mi>&amp;omega;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo> </mrow> </mrow>
那么谱表示模拟公式为：
<mrow> <msub> <mi>x</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>Re</mi> <mo>{</mo> <mn>2</mn> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>j</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msub> <mi>H</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <msqrt> <mrow> <mi>&amp;Delta;</mi> <mi>&amp;omega;</mi> </mrow> </msqrt> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>l</mi> </msub> <mi>t</mi> <mo>+</mo> <msub> <mi>&amp;phi;</mi> <mrow> <mi>k</mi> <mi>l</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </msup> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo> </mrow> </mrow>
其中，△ω＝ωu/N是频率分辨率，ωu是截止频率，N是离散频率数，ωl＝l△ω，Vkl是在区间[0,2π]均匀分布的独立随机相位角；
步骤3：将得到的三角矩阵各个时变谱作特征正交分解：
将三角分解得到的时变谱Hjk(ω,t)进行特征正交分解后，再近似表达为多个时间函数和频率函数乘积之和：
<mrow> <msub> <mi>H</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;ap;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>N</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> </munderover> <msubsup> <mi>a</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>&amp;Phi;</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>)</mo> </mrow> <mo>,</mo> <mi>j</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>n</mi> <mo>;</mo> <mi>j</mi> <mo>&amp;GreaterEqual;</mo> <mi>k</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>9</mn> <mo>)</mo> </mrow> </mrow>
其中：是频率相关矩阵Rjk的第q个特征向量；是通过计算的第q个主坐标；是包含绝大部分能量的有效项数；
步骤4：使用FFT技术模拟非平稳过程：
利用特征正交分解拟合式，将所述谱表示模拟公式，即式(8)转换为FFT直接可用的形式，如下所示：
<mrow> <msub> <mi>x</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>Re</mi> <mo>{</mo> <mn>2</mn> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>j</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>N</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> </munderover> <msubsup> <mi>a</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msubsup> <mi>&amp;Phi;</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <msqrt> <mrow> <mi>&amp;Delta;</mi> <mi>&amp;omega;</mi> </mrow> </msqrt> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>l</mi> </msub> <mi>t</mi> <mo>+</mo> <msub> <mi>&amp;phi;</mi> <mrow> <mi>k</mi> <mi>l</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </msup> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>10</mn> <mo>)</mo> </mrow> </mrow>
再采用FFT技术快速模拟。
2.根据权利要求1所述的任意非平稳随机过程快速模拟方法，其特征在于，所述步骤3中，对时变谱Hjk(ω,t)采用特征正交分解的具体过程下：
若使特征正交分解可使用，则Hjk(ω,t)需在时频域上离散化，该连续的二元函数离散为一个N维列向量：
<mrow> <msub> <mi>H</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <msup> <mrow> <mo>&amp;lsqb;</mo> <msubsup> <mi>H</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> <mn>1</mn> </msubsup> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>,</mo> <msubsup> <mi>H</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msubsup> <mi>H</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> <mi>N</mi> </msubsup> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> <mi>p</mi> <mo>=</mo> <mn>0</mn> <mo>,</mo> <mn>1</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msub> <mi>N</mi> <mi>p</mi> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>11</mn> <mo>)</mo> </mrow> </mrow>
其中，表示第l个频率点处的时间序列，通过特征正交分解，找到一组最优正交基为在这组正交基上，Hjk(p△t)的投影最大；这组基向量通过以下特征向量分解来获得：
<mrow> <msub> <mi>R</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <msubsup> <mi>&amp;Phi;</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mo>=</mo> <msubsup> <mi>&amp;lambda;</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <msubsup> <mi>&amp;Phi;</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mo>,</mo> <mi>q</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>N</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>
其中，
<mrow> <msub> <mi>R</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <msub> <mi>N</mi> <mi>p</mi> </msub> </mfrac> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mn>1</mn> </mrow> <msub> <mi>N</mi> <mi>p</mi> </msub> </munderover> <msub> <mi>H</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <msub> <mi>H</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <msup> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>13</mn> <mo>)</mo> </mrow> </mrow>
其为时间平均的频率相关矩阵；是第q个特征向量，为频率的函数，是第q个特征值；当特征值按降序排列时，低阶特征值包含更多能量，则Hjk(p△t)近似表示为：
<mrow> <msub> <mi>H</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;ap;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>N</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> </munderover> <msubsup> <mi>a</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>&amp;Phi;</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mo>,</mo> <msubsup> <mi>N</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mo>&lt;</mo> <mi>N</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>14</mn> <mo>)</mo> </mrow> </mrow>
其中，是第q个主坐标或时间函数，用确定。
3.根据权利要求1所述的任意非平稳随机过程快速模拟方法，其特征在于，所述步骤4中，采用FFT技术对谱表示模拟公式进行快速计算，其具体处理如下：
为了便于FFT应用于模拟，模拟公式改写为如下形式：
<mrow> <mtable> <mtr> <mtd> <mrow> <msub> <mi>x</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>Re</mi> <mo>{</mo> <mn>2</mn> <mi>exp</mi> <mo>&amp;lsqb;</mo> <mo>-</mo> <mi>i</mi> <mi>&amp;Delta;</mi> <mi>&amp;omega;</mi> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>&amp;rsqb;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>j</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>N</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> </munderover> <msubsup> <mi>a</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>M</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msubsup> <mi>B</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> <mi>l</mi> </mrow> </msubsup> <mo>&amp;CenterDot;</mo> <mi>exp</mi> <mo>&amp;lsqb;</mo> <mo>-</mo> <mi>i</mi> <mrow> <mo>(</mo> <mi>l</mi> <mi>&amp;Delta;</mi> <mi>&amp;omega;</mi> <mo>)</mo> </mrow> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mo>}</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>p</mi> <mo>=</mo> <mn>0</mn> <mo>,</mo> <mn>1</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>M</mi> <mo>-</mo> <mn>1</mn> </mrow> </mtd> </mtr> </mtable> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>
其中，M是时间离散点数；为：
<mrow> <msubsup> <mi>B</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> <mi>l</mi> </mrow> </msubsup> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msubsup> <mi>&amp;Phi;</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>l</mi> <mi>&amp;Delta;</mi> <mi>&amp;omega;</mi> <mo>+</mo> <mi>&amp;Delta;</mi> <mi>&amp;omega;</mi> <mo>)</mo> </mrow> <msqrt> <mrow> <mi>&amp;Delta;</mi> <mi>&amp;omega;</mi> </mrow> </msqrt> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>i&amp;phi;</mi> <mrow> <mi>k</mi> <mi>l</mi> </mrow> </msub> </mrow> </msup> </mrow> </mtd> <mtd> <mrow> <mi>l</mi> <mo>=</mo> <mn>0</mn> <mo>,</mo> <mn>1</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>N</mi> <mo>-</mo> <mn>1</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mn>0</mn> </mtd> <mtd> <mrow> <mi>l</mi> <mo>=</mo> <mi>N</mi> <mo>,</mo> <mi>N</mi> <mo>+</mo> <mn>1</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>M</mi> <mo>-</mo> <mn>1</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> </mrow>
时间和频率分辨率满足：
M△t＝2π/△ω (17)
为了避免混叠，时间步长也需要满足：
△t≤2π/(2ωu) (18)
也即
M≥2N (19)
进一步改写为：
<mrow> <msub> <mi>x</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>Re</mi> <mo>{</mo> <mn>2</mn> <mi>exp</mi> <mo>&amp;lsqb;</mo> <mo>-</mo> <mi>i</mi> <mi>&amp;Delta;</mi> <mi>&amp;omega;</mi> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>&amp;rsqb;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>j</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>N</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> </munderover> <msubsup> <mi>a</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>p</mi> <mi>&amp;Delta;</mi> <mi>t</mi> <mo>)</mo> </mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>0</mn> </mrow> <mrow> <mi>M</mi> <mo>-</mo> <mn>1</mn> </mrow> </munderover> <msubsup> <mi>B</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>k</mi> <mi>l</mi> </mrow> </msubsup> <mo>&amp;times;</mo> <mi>exp</mi> <mo>&amp;lsqb;</mo> <mo>-</mo> <mi>i</mi> <mi>l</mi> <mi>p</mi> <mfrac> <mrow> <mn>2</mn> <mi>&amp;pi;</mi> </mrow> <mi>M</mi> </mfrac> <mo>)</mo> <mo>&amp;rsqb;</mo> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow>
为充分利用FFT技术，时间的离散数满足M＝2u，u为正整数。
4.一种时不变相干的非平稳随机过程快速模拟方法，其特征在于，包括以下步骤：
步骤A：获取时不变相干的非平稳随机过程的目标功率谱和相干函数：
时不变相干的非平稳过程的演化功率谱矩阵表示为：
S(ω,t)＝D(ω,t)Γ(ω)DT(ω,t) (21)
其中
<mrow> <mi>D</mi> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>d</mi> <mi>i</mi> <mi>a</mi> <mi>g</mi> <mo>&amp;lsqb;</mo> <msqrt> <mrow> <msub> <mi>S</mi> <mn>11</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </msqrt> <mo>,</mo> <msqrt> <mrow> <msub> <mi>S</mi> <mn>22</mn> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </msqrt> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msqrt> <mrow> <msub> <mi>S</mi> <mrow> <mi>n</mi> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </msqrt> <mo>&amp;rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>22</mn> <mo>)</mo> </mrow> </mrow>
其中，ω是圆频率，γjk(ω)是随机过程xj(t)和xk(t)之间的时不变的相干函数，它们之间的相关函数表示为：
<mrow> <msub> <mi>R</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>,</mo> <mi>t</mi> <mo>+</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> <mo>=</mo> <msubsup> <mo>&amp;Integral;</mo> <mrow> <mo>-</mo> <mi>&amp;infin;</mi> </mrow> <mi>&amp;infin;</mi> </msubsup> <msqrt> <mrow> <msub> <mi>S</mi> <mrow> <mi>j</mi> <mi>j</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> <msub> <mi>S</mi> <mrow> <mi>k</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>+</mo> <mi>&amp;tau;</mi> <mo>)</mo> </mrow> </mrow> </msqrt> <msub> <mi>&amp;gamma;</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mi>i</mi> <mi>&amp;omega;</mi> <mi>&amp;tau;</mi> </mrow> </msup> <mi>d</mi> <mi>&amp;omega;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>24</mn> <mo>)</mo> </mrow> </mrow>
其中，Sjj(ω,t)为随机过程xj(t)的自谱，Skk(ω,t+t)为随机过程xk(t)的自谱；
步骤B：对相干矩阵进行三角分解，得到谱表示模拟公式：
相干矩阵Γ(ω)为正定的Hermitian矩阵，分解为：
Γ(ω)＝B(ω)BT*(ω) (25)
其中，B(ω)为下三角矩阵，表示为：
那么谱表示模拟公式为：
<mrow> <msub> <mi>x</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>Re</mi> <mo>{</mo> <mn>2</mn> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>j</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msqrt> <mrow> <msub> <mi>S</mi> <mrow> <mi>j</mi> <mi>j</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>l</mi> </msub> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </msqrt> <msub> <mi>&amp;beta;</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <msqrt> <mrow> <mi>&amp;Delta;</mi> <mi>&amp;omega;</mi> </mrow> </msqrt> <msup> <mi>e</mi> <mrow> <mo>-</mo> <mi>i</mi> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>l</mi> </msub> <mi>t</mi> <mo>+</mo> <msub> <mi>&amp;phi;</mi> <mrow> <mi>k</mi> <mi>l</mi> </mrow> </msub> <mo>)</mo> </mrow> </mrow> </msup> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>27</mn> <mo>)</mo> </mrow> </mrow>
其中，△ω＝ωu/N是频率分辨率，ωu是截止频率，N是离散频率数，ωl＝l△ω，φkl是在区间[0,2π]均匀分布的独立随机相位角；
步骤C：对关于自谱的元素作特征正交分解：
将元素近似表达为多个时间函数和频率函数乘积之和：
<mrow> <msqrt> <mrow> <msub> <mi>S</mi> <mrow> <mi>j</mi> <mi>j</mi> </mrow> </msub> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>,</mo> <mi>t</mi> <mo>)</mo> </mrow> </mrow> </msqrt> <mo>&amp;ap;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>N</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>j</mi> </mrow> </msubsup> </munderover> <msubsup> <mi>a</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>j</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <msubsup> <mi>&amp;Phi;</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>j</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>28</mn> <mo>)</mo> </mrow> </mrow>
其中：是频率相关矩阵Rjj的第q个特征向量；是通过计算的第q个主坐标；是包含绝大部分能量的有效项数；
步骤D：优化模拟公式，并采用FFT技术模拟非平稳过程：
利用特征正交分解拟合式，将所述谱表示模拟公式，即式(27)转换为FFT直接可用的形式，然后交换公式求和的顺序，减少FFT的执行次数，如下所示：
<mrow> <msub> <mi>x</mi> <mi>j</mi> </msub> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <mo>=</mo> <mi>Re</mi> <mo>{</mo> <mn>2</mn> <msqrt> <mrow> <mi>&amp;Delta;</mi> <mi>&amp;omega;</mi> </mrow> </msqrt> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mn>1</mn> </mrow> <msubsup> <mi>N</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>j</mi> </mrow> </msubsup> </munderover> <msubsup> <mi>a</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>j</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>)</mo> </mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <msubsup> <mi>&amp;Phi;</mi> <mi>q</mi> <mrow> <mi>j</mi> <mi>j</mi> </mrow> </msubsup> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>i&amp;omega;</mi> <mi>l</mi> </msub> <mi>t</mi> </mrow> </msup> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>k</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>j</mi> </munderover> <msub> <mi>&amp;beta;</mi> <mrow> <mi>j</mi> <mi>k</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>&amp;omega;</mi> <mi>l</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>e</mi> <mrow> <mo>-</mo> <msub> <mi>i&amp;phi;</mi> <mrow> <mi>k</mi> <mi>l</mi> </mrow> </msub> </mrow> </msup> <mo>}</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>29</mn> <mo>)</mo> </mrow> </mrow>
其中，△ω＝ωu/N是频率分辨率，ωu是截止频率，N是离散频率数，ωl＝l△ω；
再采用FFT技术快速模拟非平稳过程时程样本。
CN201710351315.4A 2017-05-18 2017-05-18 一种基于特征正交分解的非平稳随机过程快速模拟方法 Expired - Fee Related CN107247687B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710351315.4A CN107247687B (zh) 2017-05-18 2017-05-18 一种基于特征正交分解的非平稳随机过程快速模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710351315.4A CN107247687B (zh) 2017-05-18 2017-05-18 一种基于特征正交分解的非平稳随机过程快速模拟方法

Publications (2)

Publication Number Publication Date
CN107247687A true CN107247687A (zh) 2017-10-13
CN107247687B CN107247687B (zh) 2021-06-08

Family

ID=60016645

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710351315.4A Expired - Fee Related CN107247687B (zh) 2017-05-18 2017-05-18 一种基于特征正交分解的非平稳随机过程快速模拟方法

Country Status (1)

Country Link
CN (1) CN107247687B (zh)

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110196960A (zh) * 2019-04-24 2019-09-03 重庆大学 一种基于二维fft的各态历经风场高效模拟方法
CN110276087A (zh) * 2019-03-12 2019-09-24 南昌大学 一种非平稳随机扭转振动疲劳寿命预测方法
CN110826197A (zh) * 2019-10-21 2020-02-21 西南交通大学 一种基于改进Cholesky分解闭合解的风速场模拟方法
CN111368392A (zh) * 2019-12-31 2020-07-03 重庆大学 一种基于memd与srm的单样本非平稳风速模拟方法
CN113806845A (zh) * 2021-09-16 2021-12-17 四川农业大学 一种基于pod插值的非平稳风场高效模拟方法
CN115510381A (zh) * 2022-09-27 2022-12-23 中国海洋大学 一种海上风机多元相干效应风场载荷构建方法
CN115526125A (zh) * 2022-09-22 2022-12-27 四川农业大学 一种基于数值截断的随机风速场高效模拟方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104408295A (zh) * 2014-11-10 2015-03-11 浙江大学 一种大跨桥梁下部结构风-浪耦合作用荷载数值模拟方法
CN104516771A (zh) * 2015-01-22 2015-04-15 黄国庆 一种非平稳随机过程高效模拟方法
CN105426594A (zh) * 2015-11-06 2016-03-23 西南交通大学 一种基于时空场和条件插值的平稳均质风场快速模拟方法
US20160171134A1 (en) * 2014-12-11 2016-06-16 Siemens Aktiengesellschaft Automatic ranking of design parameter significance for fast and accurate cae-based design space exploration using parameter sensitivity feedback
CN106682277A (zh) * 2016-12-06 2017-05-17 西南交通大学 一种非平稳随机过程快速模拟方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104408295A (zh) * 2014-11-10 2015-03-11 浙江大学 一种大跨桥梁下部结构风-浪耦合作用荷载数值模拟方法
US20160171134A1 (en) * 2014-12-11 2016-06-16 Siemens Aktiengesellschaft Automatic ranking of design parameter significance for fast and accurate cae-based design space exploration using parameter sensitivity feedback
CN104516771A (zh) * 2015-01-22 2015-04-15 黄国庆 一种非平稳随机过程高效模拟方法
CN105426594A (zh) * 2015-11-06 2016-03-23 西南交通大学 一种基于时空场和条件插值的平稳均质风场快速模拟方法
CN106682277A (zh) * 2016-12-06 2017-05-17 西南交通大学 一种非平稳随机过程快速模拟方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
GUOQING HUANG: "Application of Proper Orthogonal Decomposition in Fast Fourier Transform—Assisted Multivariate Nonstationary Process Simulation", 《JOURNAL OF ENGINEERING MECHANICS》 *
GUOQING HUANG等: "New formulation of Cholesky decomposition and applications in stochastic simulation", 《PROBABILISTIC ENGINEERING MECHANICS》 *
李锦华等: "《下击暴流非平稳脉动风速数值模拟》", 《振动与冲击》 *

Cited By (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110276087A (zh) * 2019-03-12 2019-09-24 南昌大学 一种非平稳随机扭转振动疲劳寿命预测方法
CN110196960A (zh) * 2019-04-24 2019-09-03 重庆大学 一种基于二维fft的各态历经风场高效模拟方法
CN110196960B (zh) * 2019-04-24 2022-12-27 重庆大学 一种基于二维fft的各态历经风场高效模拟方法
CN110826197A (zh) * 2019-10-21 2020-02-21 西南交通大学 一种基于改进Cholesky分解闭合解的风速场模拟方法
CN110826197B (zh) * 2019-10-21 2022-11-01 西南交通大学 一种基于改进Cholesky分解闭合解的风速场模拟方法
CN111368392A (zh) * 2019-12-31 2020-07-03 重庆大学 一种基于memd与srm的单样本非平稳风速模拟方法
CN111368392B (zh) * 2019-12-31 2024-04-05 重庆大学 一种基于memd与srm的单样本非平稳风速模拟方法
CN113806845A (zh) * 2021-09-16 2021-12-17 四川农业大学 一种基于pod插值的非平稳风场高效模拟方法
CN115526125A (zh) * 2022-09-22 2022-12-27 四川农业大学 一种基于数值截断的随机风速场高效模拟方法
CN115526125B (zh) * 2022-09-22 2023-09-19 四川农业大学 一种基于数值截断的随机风速场高效模拟方法
CN115510381A (zh) * 2022-09-27 2022-12-23 中国海洋大学 一种海上风机多元相干效应风场载荷构建方法
CN115510381B (zh) * 2022-09-27 2023-08-22 中国海洋大学 一种海上风机多元相干效应风场载荷构建方法

Also Published As

Publication number Publication date
CN107247687B (zh) 2021-06-08

Similar Documents

Publication Publication Date Title
CN107247687A (zh) 一种基于特征正交分解的非平稳随机过程快速模拟方法
Liu et al. Dimension reduction of Karhunen-Loeve expansion for simulation of stochastic processes
Heitmann et al. The coyote universe. I. Precision determination of the nonlinear matter power spectrum
Chang Estimation of wind energy potential using different probability density functions
CN103139907A (zh) 一种利用指纹法的室内无线定位方法
CN106202977B (zh) 一种基于盲源分离算法的低频振荡模式分析方法
CN110083920B (zh) 一种地震作用下非比例阻尼结构随机响应的分析方法
Liu et al. Validity analysis of maximum entropy distribution based on different moment constraints for wind energy assessment
CN104408295B (zh) 一种大跨桥梁下部结构风‑浪耦合作用荷载数值模拟方法
Bao et al. Fast simulation of non-stationary wind velocity based on time-frequency interpolation
CN106202734A (zh) 基于高斯径向基函数的全局灵敏度分析方法
Le et al. Computationally efficient stochastic approach for the fragility analysis of vertical structures subjected to thunderstorm downburst winds
CN113806845A (zh) 一种基于pod插值的非平稳风场高效模拟方法
Dai et al. Difference-based variance estimation in nonparametric regression with repeated measurement data
Zlatev et al. Richardson Extrapolation combined with the sequential splitting procedure and the θ-method
CN104239598A (zh) 一种面向动态系统模型验证的多元数据分析方法
CN101364245B (zh) 多极子数据库的电磁环境预测系统
CN113609700A (zh) 一种山区桥梁完全非平稳风场的模拟方法
CN109033181B (zh) 一种复杂地形地区风场地理数值模拟方法
CN104036118B (zh) 一种电力系统并行化轨迹灵敏度获取方法
García Vera et al. Multilevel algorithm for flow observables in gauge theories
CN105960613A (zh) 系统辨识装置
Yoon et al. On the propagation of information and the use of localization in ensemble Kalman filtering
JP7156613B2 (ja) 津波予測装置、方法、及びプログラム
CN112131638A (zh) 大跨屋盖结构的风致动力特性类型判定方法及终端设备

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210608