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

本发明涉及随机信号模拟
技术领域：
，具体为一种基于特征正交分解的非平稳随机过程快速模拟方法。
背景技术：
：地震和瞬时极端风等毁灭性的激励表现出非平稳特性，通常用非平稳过程来描述。当涉及到非线性、系统随机性、参数激励和其他某些随机问题时，非平稳过程的蒙特卡洛模拟(mcs)是结构响应分析中的一项重要任务。此外，mcs是评价其他随机方法准确性的基准。基于演化功率谱密度函数(epsdevolutionarypowerspectraldensity)的谱表示方法，由于其准确性和简便性被广泛应用于随机过程的模拟中。然而，epsd的时变特性使模拟出现了两个困难。一方面，功率谱矩阵在每个时刻和频率都要进行一次三角分解。另一方面，由于不能使用快速傅里叶变换(fft)，非平稳随机过程的模拟效率较为低下。(li,y.,andkareem,a.simulationofmultivariatenonstationaryrandomprocessesbyfft[j].journalofengineeringmechanics,1991,117,1037-58.)使用了三角函数和勒让德多项式函数来拟合分解后的epsd矩阵，建立了使用fft的谱表示方法。(huang,g.anefficientsimulationapproachformultivariatenonstationaryprocess:hybridofwaveletandspectralrepresentationmethod[j].probabilistic.engineeringmechanics.,2014,37,74-83.)使用小波解耦时变epsd矩阵，发展了一种基于快速傅里叶变换的谱表示方法。(huang,g.applicationofproperorthogonaldecompositioninfastfouriertransform—assistedmultivariatenonstationaryprocesssimulation[j].journalofengineeringmechanics2015,141.7:04015015.)采用特征正交分解(podproperorthogonaldecomposition)来分离分解后的epsd矩阵，并利用快速傅里叶变换来提高模拟速度。这里的基于pod的分解是基于数据的，最优的并且是易于工程应用的。(peng,l.,huang,g.,chen,x.,andkareem,a.simulationofmultivariatenonstationaryrandomprocess:ahybridstochasticwaveandproperorthogonaldecompositionapproach[j].journalofengineeringmechanics(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.applicationofproperorthogonaldecompositioninfastfouriertransform—assistedmultivariatenonstationaryprocesssimulation[j].journalofengineeringmechanics141.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)andh55(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分解32611736003600pod15926736442fft7578总时间19251241021645此外，通过进一步计算，当有50、175、301个模拟点时，总的时间比分别是117、388、645。可以发现，模拟点数越多，效率提升越明显。当前第1页12