本發(fā)明涉及隨機信號模擬
技術領域：
，具體為一種基于特征正交分解的非平穩(wěn)隨機過程快速模擬方法。
背景技術：
：地震和瞬時極端風等毀滅性的激勵表現(xiàn)出非平穩(wěn)特性，通常用非平穩(wěn)過程來描述。當涉及到非線性、系統(tǒng)隨機性、參數(shù)激勵和其他某些隨機問題時，非平穩(wěn)過程的蒙特卡洛模擬(mcs)是結構響應分析中的一項重要任務。此外，mcs是評價其他隨機方法準確性的基準?；谘莼β首V密度函數(shù)(epsdevolutionarypowerspectraldensity)的譜表示方法，由于其準確性和簡便性被廣泛應用于隨機過程的模擬中。然而，epsd的時變特性使模擬出現(xiàn)了兩個困難。一方面，功率譜矩陣在每個時刻和頻率都要進行一次三角分解。另一方面，由于不能使用快速傅里葉變換(fft)，非平穩(wěn)隨機過程的模擬效率較為低下。(li,y.,andkareem,a.simulationofmultivariatenonstationaryrandomprocessesbyfft[j].journalofengineeringmechanics,1991,117,1037-58.)使用了三角函數(shù)和勒讓德多項式函數(shù)來擬合分解后的epsd矩陣，建立了使用fft的譜表示方法。(huang,g.anefficientsimulationapproachformultivariatenonstationaryprocess:hybridofwaveletandspectralrepresentationmethod[j].probabilistic.engineeringmechanics.,2014,37,74-83.)使用小波解耦時變epsd矩陣，發(fā)展了一種基于快速傅里葉變換的譜表示方法。(huang,g.applicationofproperorthogonaldecompositioninfastfouriertransform—assistedmultivariatenonstationaryprocesssimulation[j].journalofengineeringmechanics2015,141.7:04015015.)采用特征正交分解(podproperorthogonaldecomposition)來分離分解后的epsd矩陣，并利用快速傅里葉變換來提高模擬速度。這里的基于pod的分解是基于數(shù)據(jù)的，最優(yōu)的并且是易于工程應用的。(peng,l.,huang,g.,chen,x.,andkareem,a.simulationofmultivariatenonstationaryrandomprocess:ahybridstochasticwaveandproperorthogonaldecompositionapproach[j].journalofengineeringmechanics(2017accepted))提出了一種新的基于隨機波和pod的混合方法，用于模擬沿直線分布的多點非平穩(wěn)過程。其中二維fft的使用大大地提高了模擬的效率?？偟膩碚f，在非平穩(wěn)模擬方面，雖然已經取得了很大的進步，但是仍然有必要提高模擬的效率和適用性，特別是對于有許多模擬點的現(xiàn)代龐大的結構。此外，一些非平穩(wěn)激勵的相干函數(shù)具有明顯的時變特征，如下?lián)舯┝黠L。因此，有必要把常規(guī)的模擬方法拓展到時變相干的非平穩(wěn)過程。技術實現(xiàn)要素：針對上述問題，本發(fā)明的目的在于提供一種可有效地解決目前各種非平穩(wěn)模擬方法存在的應用局限以及效率低下等問題的基于特征正交分解的非平穩(wěn)隨機過程快速模擬方法，采用該方法模擬出的樣本準確，模擬效率高的。技術方案如下：一種任意非平穩(wěn)隨機過程快速模擬方法，包括以下步驟：步驟1：獲取時變相干的多點非平穩(wěn)隨機過程的目標功率譜和相干函數(shù)：獲取零均值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)的時變的相干函數(shù)，由相干函數(shù)組成的相干矩陣γ(ω,t)由下式給出：步驟2：對演化功率譜矩陣進行三角分解，得到分解時變譜，進而得到譜表示模擬公式：通過cholesky分解，演化功率譜矩陣分解為：其中，h(ω,t)為下三角矩陣形，表示為：任意兩個隨機過程之間的演化功率譜sjk(ω,t)表達為：它們的相關函數(shù)表達為：那么譜表示模擬公式為：其中，△ω＝ωu/n是頻率分辨率，ωu是截止頻率，n是離散頻率數(shù)，ωl＝l△ω，φkl是在區(qū)間[0,2π]均勻分布的獨立隨機相位角；步驟3：將得到的三角矩陣各個時變譜作特征正交分解：將三角分解得到的時變譜hjk(ω,t)進行特征正交分解后，再近似表達為多個時間函數(shù)和頻率函數(shù)乘積之和：其中：是頻率相關矩陣rjk的第q個特征向量；是通過計算的第q個主坐標；是包含絕大部分能量的有效項數(shù)；步驟4：使用fft技術模擬非平穩(wěn)過程：利用特征正交分解擬合式，將所述譜表示模擬公式，即式(8)轉換為fft直接可用的形式，如下所示：再采用fft技術快速模擬。進一步的，所述步驟3中，對時變譜hjk(ω,t)采用特征正交分解的具體過程下：若使特征正交分解可使用，則hjk(ω,t)需在時頻域上離散化，該連續(xù)的二元函數(shù)離散為一個n維列向量：其中，表示第l個頻率點處的時間序列，通過特征正交分解，找到一組最優(yōu)正交基為在這組正交基上，hjk(p△t)的投影最大；這組基向量通過以下特征向量分解來獲得：其中，其為時間平均的頻率相關矩陣；是第q個特征向量，為頻率的函數(shù)，是第q個特征值；當特征值按降序排列時，低階特征值包含更多能量，則hjk(p△t)近似表示為：其中，是第q個主坐標或時間函數(shù)，用確定。更進一步的，所述步驟4中，采用fft技術對譜表示模擬公式進行快速計算，其具體處理如下：為了便于fft應用于模擬，模擬公式改寫為如下形式：其中，m是時間離散點數(shù)；為：時間和頻率分辨率滿足：m△t＝2π/△ω(17)為了避免混疊，時間步長也需要滿足：△t≤2π/(2ωu)(18)也即m≥2n(19)進一步改寫為：為充分利用fft技術，時間的離散數(shù)滿足m＝2u，u為正整數(shù)。一種時不變相干的非平穩(wěn)隨機過程快速模擬方法，包括以下步驟：步驟a：獲取時不變相干的非平穩(wěn)隨機過程的目標功率譜和相干函數(shù)：時不變相干的非平穩(wěn)過程的演化功率譜矩陣表示為：s(ω,t)＝d(ω,t)γ(ω)dt(ω,t)(21)其中其中，ω是圓頻率，γjk(ω)是隨機過程xj(t)和xk(t)之間的時不變的相干函數(shù)，它們之間的相關函數(shù)表示為：其中，sjj(ω,t)為隨機過程xj(t)的自譜，skk(ω,t+t)為隨機過程xk(t)的自譜；步驟b：對相干矩陣進行三角分解，得到譜表示模擬公式：相干矩陣γ(ω)為正定的hermitian矩陣，分解為：其中，b(ω)為下三角矩陣，表示為：那么譜表示模擬公式為：其中，△ω＝ωu/n是頻率分辨率，ωu是截止頻率，n是離散頻率數(shù)，ωl＝l△ω，φkl是在區(qū)間[0,2π]均勻分布的獨立隨機相位角；步驟c：對關于自譜的元素作特征正交分解，得到譜表示模擬公式：將元素近似表達為多個時間函數(shù)和頻率函數(shù)乘積之和：其中：是頻率相關矩陣rjj的第q個特征向量；是通過計算的第q個主坐標；是包含絕大部分能量的有效項數(shù)；步驟d：優(yōu)化模擬公式，并采用fft技術模擬非平穩(wěn)過程：利用特征正交分解擬合式，將所述譜表示模擬公式，即式(27)轉換為fft直接可用的形式，然后交換公式求和的順序，減少fft的執(zhí)行次數(shù)，如下所示：其中，△ω＝ωu/n是頻率分辨率，ωu是截止頻率，n是離散頻率數(shù)，ωl＝l△ω；再采用fft技術快速模擬非平穩(wěn)過程時程樣本。本發(fā)明的有益效果是：本發(fā)明的方法能夠模擬任意的非平穩(wěn)過程，由于pod擬合的精度高，所以模擬出的樣本準確；三角分解及fft技術的優(yōu)化使用使模擬效率大大提升；有效地解決了目前各種非平穩(wěn)模擬方法存在的應用局限以及效率低下等問題，具有廣闊的應用前景。附圖說明圖1a為實例1中的模擬點3和模擬點5分解后的時變譜的主坐標圖1b為實例1中的模擬點3和模擬點5分解后的特征向量。圖2a為實例1中模擬點3采用兩種模擬方法擬合的時變譜與目標值的對比。圖2b為實例1中模擬點5采用兩種模擬方法擬合的時變譜與目標值的對比。圖3為實例1中各模擬點的樣本時程。圖4為實例1中模擬點3自相關函數(shù)與目標值的對比。圖5為實例1中模擬點5的自相關函數(shù)與目標值的對比。圖6為實例1中模擬點3模擬點5的互相關函數(shù)與目標值的對比。具體實施方式下面結合附圖和具體實施方式對本發(fā)明作進一步詳細的說明。本實施例示出一種基于特征正交分解(pod)采用快速傅立葉變換(fft)的譜表示模擬方法：對于空間任意分布的、時變相干的多點非平穩(wěn)過程，首先對演化功率譜矩陣進行三角分解，然后使用特征正交分解(pod)將三角矩陣的各個時頻耦合的元素分別表示為若干個時間函數(shù)和頻率函數(shù)乘積之和，最后采用fft技術進行快速模擬。特別地，當該方法應用于空間任意分布的、時不變相干的多點非平穩(wěn)過程時，首先對相干矩陣進行三角分解，然后使用pod將自譜組成的對角矩陣的時頻耦合的元素，分別表示為若干個時間函數(shù)和頻率函數(shù)乘積之和，最后優(yōu)化模擬公式，用極少次的fft技術進行高效模擬。一種任意的非平穩(wěn)隨機過程快速模擬方法，具體步驟如下：步驟1：獲取多點非平穩(wěn)隨機過程的目標功率譜和相干函數(shù)：獲取零均值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)的時變的相干函數(shù)。由相干函數(shù)組成的相干矩陣γ(ω,t)由下式給出：然后，對頻率和時間離散化處理，進而獲得在離散時頻域上的演化功率譜矩陣。步驟2：在離散的時頻點上對演化功率譜矩陣進行三角分解，得到離散的分解時變譜。通過cholesky分解，epsd矩陣可以作出以下分解：其中：h(ω,t)是一個下三角矩陣形如：它的對角元素一般是非負實數(shù)；而非對角元素通常是復數(shù)。任意兩個隨機過程之間的演化功率譜sjk(ω,t)表達為：它們的相關函數(shù)表達為：那么譜表示模擬公式為：其中，△ω＝ωu/n是頻率分辨率，ωu是截止頻率，n是離散頻率數(shù)，ωl＝l△ω，φkl是在區(qū)間[0,2π]均勻分布的獨立隨機相位角。步驟3：將得到的三角矩陣各個離散的時變譜作特征正交分解：將分解后的時變譜采用特征正交分解(pod)的具體處理如下：為了使pod能夠使用，hjk(ω,t)應該在時頻域上離散化。該連續(xù)的二元函數(shù)可以離散為一個n維列向量：其中表示第l個頻率點處的時間序列。通過pod，可以找到一組最優(yōu)正交基為在這組正交基上，hjk(p△t)的投影是最大的。這組基向量可以通過以下特征向量分解來獲得，其中它可以被解釋為時間平均的頻率相關矩陣；是第q個特征向量，它也是一個頻率的函數(shù)，是第q個特征值。當特征值按降序排列時，低階特征值包含更多能量。那么，hjk(p△t)可以近似表示為：其中，是第q個主坐標(或時間函數(shù))，可以用確定。將分解得到的時變譜hjk(ω,t)進一步近似表達為多個時間序列和頻率序列乘積之和，如下所示：其中是頻率相關矩陣rjk的第q個特征向量；是通過計算的第q個主坐標；是選取的包含絕大部分能量的有效項數(shù)。利用pod擬合式(42)，將譜表示模擬公式(37)轉換為fft直接可用的形式：步驟4：使用fft技術高效地模擬非平穩(wěn)過程：式(43)可進一步寫為離散化的形式：其中m是時間離散點數(shù)；△t是時間分辨率，時間分辨率△t和頻率分辨率△ω滿足：m△t＝2π/△ω(45)為了避免混疊，時間步長也需要滿足：△t≤2π/(2ωu)(46)也即m≥2n(47)進一步改寫為：其中，滿足m＝2u，u是一個正整數(shù)；為：其中△ω＝ωu/n是頻率分辨率，ωu是截止頻率，n是離散頻率數(shù)。然后使用fft技術快速模擬出非平穩(wěn)時程樣本。進一步，通過下式確定目標相關函數(shù)：然后，將模擬的樣本計算的相關函數(shù)與目標相關函數(shù)進行對比，驗證模擬的精度。值得指出的是，經典的譜表示模擬公式將分解后的時變譜表示為極坐標的形式，如下：其中|·|表示復數(shù)的模；θjk(ω,t)是hjk(ω,t)的相位角，通過下式確定：θjk(ω,t)＝tan-1{im[hjk(ω,t)]/re[hjk(ω,t)]}(52)其中：im和re分別表示虛部和實部。但是這樣求得的相位角只在(-π/2,π/2)的角度范圍，而實際的相位角應在(-π,π)的范圍內。但是，本發(fā)明所提方法避免了這樣的問題。對于時不變相干的非平穩(wěn)過程的模擬，其具體實施步驟為：步驟a：獲取時不變非平穩(wěn)隨機過程的目標功率譜和相干函數(shù)：時不變相干的非平穩(wěn)過程的演化功率譜矩陣可以表示為：s(ω,t)＝d(ω,t)γ(ω)dt(ω,t)(53)其中其中：γjk(ω)是xj(t)和xk(t)之間的時不變的相干函數(shù)。然后，對頻率和時間離散化處理，進而獲得在離散時頻域上的功率譜和相干函數(shù)。它們之間的相關函數(shù)表示為：其中，sjj(ω,t)為隨機過程xj(t)的自譜，skk(ω,t+t)為隨機過程xk(t)的自譜。步驟b：對離散頻域上的相干矩陣進行三角分解：相干矩陣γ(ω)也是一個正定的hermitian矩陣，它可以被分解為其中b(ω)是一個下三角矩陣，可以表示為那么譜表示模擬公式為：其中，△ω＝ωu/n是頻率分辨率，ωu是截止頻率，n是離散頻率數(shù)，ωl＝l△ω，φkl是在區(qū)間[0,2π]均勻分布的獨立隨機相位角。步驟c：對關于自譜的元素作特征正交分解：將元素近似表達為多個時間序列和頻率序列乘積之和，如下所示：其中：是頻率相關矩陣rjj的第q個特征向量；是通過計算的第q個主坐標；是包含絕大部分能量的有效項數(shù)。利用上述pod擬合式(60)，將譜表示模擬公式(59)轉換為fft直接可用的形式：步驟d：優(yōu)化模擬公式，并使用fft技術高效地模擬非平穩(wěn)過程：模擬公式(61)可進一步寫為離散化的形式，然后交換公式求和的順序，減少fft的執(zhí)行次數(shù)，如下所示：其中，最后使用fft技術快速模擬非平穩(wěn)過程的時程樣本。進一步目標相關函數(shù)可確定為：然后，將模擬的樣本計算的相關函數(shù)與目標相關函數(shù)進行對比，驗證模擬的精度。下面通過具體的實例對本發(fā)明以及本發(fā)明效果進行進一步的說明。2015年,(huang,guoqing.applicationofproperorthogonaldecompositioninfastfouriertransform—assistedmultivariatenonstationaryprocesssimulation[j].journalofengineeringmechanics141.7(2015):04015015.)將pod應用于非平穩(wěn)隨機過程的譜表示方法，提出了一種可以采用fft的非平穩(wěn)隨機過程快速模擬方法。以下通過兩個實例分別與huang的方法進行對比，驗證本發(fā)明方法的精度和效率。實例1：模擬精度驗證本例將討論沿高層建筑下?lián)舯┝黠L場的模擬。首先，通過在0.9、2.4、4、10、116、158和200m高度處(編號1-7)實測的下?lián)舯┝黠L速數(shù)據(jù)，估計出該多點非平穩(wěn)過程的演化功率譜和相干函數(shù)。從相干函數(shù)的特性發(fā)現(xiàn)它是時變的，然后按照以上時變相干的非平穩(wěn)過程的模擬步驟進行樣本生成。圖1a和圖1b分別顯示了點3和點5對應的分解后的時變譜h33(f,t)andh55(f,t)的前兩階主坐標和特征向量。可以發(fā)現(xiàn)：相應于較小特征值的主坐標(時間函數(shù))包含了較大的能量。圖2a和圖2b比較了所提方法擬合的分解后的譜與相應的目標值。其中，使用pod時，僅僅采用前5階就保留了99％的能量。因此，該方法在匹配分解后的譜方面具有優(yōu)秀的表現(xiàn)。為了便于比較，由huang的方法擬合的分解后的譜也繪制于圖2a和圖2b中，這種方法使用了前32階進行pod擬合。從圖中可以看到，這種方法在擬合分解后的譜時可能會帶來較大誤差。在模擬中，截止頻率被取為0.5hz、頻率增量為0.00049hz。生成的時程樣本的時間分辨率為1s，持續(xù)時間為2048s。圖3顯示了模擬的7個點的時程樣本。然后用10000個模擬樣本計算自相關/互相關函數(shù)。圖4、圖5、圖6對比了在0、4、8秒三個時間差上，估算的相關函數(shù)與目標值?？梢钥闯?，本發(fā)明具有很高的模擬精度。實例2：模擬效率對比一座主跨為2000米，邊跨為1000米的大跨度懸索橋沿橋面的風場被研究，為了簡單起見，只考慮垂直于橋軸線方向的風速分量。邊跨和主跨上兩個相鄰的模擬點之間的距離分別取為10米和20米。因此，共301個模擬點的風場將會被模擬。該風場特性根據(jù)實測的臺風記錄構造。同前面的例子一樣，從實測的數(shù)據(jù)估計出臺風的功率譜。然后，假定橋梁中跨模擬點的自譜等于該臺風譜，而邊跨各點的自譜取為0.8倍的臺風譜。此外，相干函數(shù)采用了davenport的指數(shù)函數(shù)模型：其中：λ＝7是衰減因子，ξ是橋面上兩點間距，u0＝40m/s是橋面高度處的平均風速。在模擬中，截止頻率取為0.5hz，頻率增量為0.000244hz。生成的時程樣本的時間分辨率為1s、持續(xù)時間為4096s。當要確保99％的能量時，pod擬合需要的近似項的數(shù)量為基于本算例我們對比了它和黃的方法的在這三方面的計算時間，如表1所示。其中，黃的方法中pod的計算時間包括了求解線性聯(lián)立方程所用的時間。表中時間比是指黃的方法的計算時間與所提方法的時間的比值。為了不失一般性，對各點自譜都使用了pod。從表中可以發(fā)現(xiàn)，總的模擬效率提高了645倍。特別地，如果只對兩個不同自譜應用pod，那么總的效率可提升3600倍。表1實例2中兩種方法模擬效率的對比(s)計算過程本發(fā)明方法huang方法時間比cholesky分解32611736003600pod15926736442fft7578總時間19251241021645此外，通過進一步計算，當有50、175、301個模擬點時，總的時間比分別是117、388、645?？梢园l(fā)現(xiàn)，模擬點數(shù)越多，效率提升越明顯。當前第1頁12