一種指定功率譜密度與峭度值的對(duì)稱(chēng)威布爾分布道路譜重構(gòu)方法

叢 楠， 陳俊達(dá)， 任焱晞， 王彬星， 李青霞

(1. 工程裝備系統(tǒng)工程研究所，北京 100093； 2. 71811部隊(duì)，湖北 孝感 432900)

功率譜密度作為路面不平度在(空間)頻域上能量分布的指征，自20世紀(jì)70年代起，已成為描述道路不平度的最主要形式[1]。當(dāng)前，無(wú)論是普通公路還是各種專(zhuān)用試驗(yàn)道路，功率譜密度均是道路不平度的最主要數(shù)值特征。在進(jìn)行道路模擬時(shí)，定義在(空間)頻域上的功率譜密度須“重構(gòu)”為隨機(jī)時(shí)域或空間域歷程，方能被輸入振動(dòng)臺(tái)架開(kāi)展試驗(yàn)。目前成熟的以諧波疊加法和線(xiàn)性濾波為代表的路譜重構(gòu)方法均僅能產(chǎn)生高斯分布的隨機(jī)歷程[2-5]。實(shí)測(cè)數(shù)據(jù)表明，高斯路形譜僅適用于路況良好的等級(jí)道路，對(duì)于路況較差或?qū)Ｓ每煽啃栽囼?yàn)道路，由于其往往包含有超過(guò)“標(biāo)準(zhǔn)”預(yù)期的不規(guī)則性，使得其幅值分布具有“厚尾”特征而呈現(xiàn)出明顯的非高斯分布特性[6-7]。

對(duì)于指定功率譜的非高斯隨機(jī)過(guò)程的模擬，當(dāng)前主要有零記憶非線(xiàn)性變換法(Zero Memory Non-linearity，ZMNL)和球不變隨機(jī)過(guò)程法(Spherically Invariant Random Process，SIRP)兩種成熟方法[8-10]。此外，蔣瑜等[11]通過(guò)在進(jìn)行IFFT時(shí)增加特定的相位信息，提出了一種生成具有指定功率譜、偏斜度和峭度的非高斯信號(hào)的方法。然而上述方法中前兩者一般只能用于生成非高斯幅度信號(hào)(即正實(shí)數(shù)信號(hào)或復(fù)數(shù)信號(hào)，而不是路譜重構(gòu)所需的正負(fù)交替且幅值概率密度呈對(duì)稱(chēng)分布的非高斯信號(hào))，后者雖然可以產(chǎn)生對(duì)稱(chēng)非高斯信號(hào)，但無(wú)法得到確定的非高斯幅值概率密度分布函數(shù)。

由于高斯隨機(jī)過(guò)程擁有固定的峭度值3而一般非高斯隨機(jī)過(guò)程其峭度值往往>3(此時(shí)稱(chēng)為超高斯隨機(jī)過(guò)程),因此，在工程實(shí)踐中，峭度值往往可作為鑒別高斯與非高斯特性的關(guān)鍵依據(jù)。本文通過(guò)推導(dǎo)威布爾分布隨機(jī)過(guò)程的分布參數(shù)與功率譜密度及峭度值的關(guān)系，并在ZMNL方法的基礎(chǔ)上，提出了一種產(chǎn)生具有指定峭度和功率譜密度，且幅度滿(mǎn)足威布爾分布的隨機(jī)信號(hào)的方法。通過(guò)將該方法應(yīng)用于試驗(yàn)場(chǎng)道路譜重構(gòu)，實(shí)現(xiàn)了在保持功率譜密度(函數(shù))曲線(xiàn)的前提下，生成具有指定峭度值的路形譜，從而為道路模擬試驗(yàn)以及各種車(chē)輛仿真試驗(yàn)提供了一種新的道路模型建模方法。

1 指定峭度的對(duì)稱(chēng)幅度威布爾分布隨機(jī)過(guò)程

威布爾分布是一類(lèi)包含了從指數(shù)分布到正態(tài)分布的分布族，其概率密度函數(shù)具有很高的包容性，因此往往能在很寬的范圍內(nèi)與試驗(yàn)數(shù)據(jù)進(jìn)行匹配。當(dāng)前，在信號(hào)模擬方面，幅度呈威布爾分布的信號(hào)已被成功應(yīng)用于模擬地物、海洋、云雨等自然條件在近距離對(duì)雷達(dá)等電子設(shè)備所造成的嚴(yán)重干擾[12]。對(duì)于路面起伏的幅度，也往往能夠用威布爾分布進(jìn)行準(zhǔn)確的描述。

對(duì)于滿(mǎn)足雙參數(shù)威布爾分布隨機(jī)過(guò)程X～Wweibull(p,q)，其具有如式(1)所示的概率密度函數(shù)：

(1)

式中：p>0為形狀參數(shù)，表明了分布函數(shù)曲線(xiàn)的偏斜情況，當(dāng)p=1時(shí)，威布爾分布等同于指數(shù)分布，當(dāng)p=2時(shí)為瑞利分布，當(dāng)p<3后則接近正態(tài)分布；q>0為尺度參數(shù)，表明了分布的中位數(shù)。

根據(jù)威布爾分布性質(zhì)，其k階原點(diǎn)矩為[13]：

(2)

式中：Γ(·)為伽馬函數(shù)。

根據(jù)峭度的定義及式(2)可知，其峭度為：

(3)

由于滿(mǎn)足上述雙參數(shù)威布爾分布的隨機(jī)變量是非負(fù)的，而實(shí)際道路起伏是有正有負(fù)，不妨令|Y|=X～Wweibull(p,q)，并令幅值概率密度函數(shù)沿0值縱軸對(duì)稱(chēng)(以下稱(chēng)該分布為對(duì)稱(chēng)幅度威布爾分布)，則其幅值的概率密度函數(shù)為：

(4)

由式(4)不難得出，由于Y的幅值概率密度的對(duì)稱(chēng)性，因此其奇數(shù)階矩(如均值和偏斜度)均為0，而偶數(shù)階原點(diǎn)矩仍與X相同。因此，式(3)退化為：

(5)

由式(5)可知，隨機(jī)過(guò)程Y的峭度僅與形狀參數(shù)p有關(guān)。當(dāng)p在區(qū)間[0.5,1.5]范圍內(nèi)取值時(shí)，峭度值的變化如圖1所示。然而，根據(jù)伽馬函數(shù)定義，當(dāng)p≠0.5, 1, 2時(shí)，由式(5)無(wú)法給出峭度K的精確解。同理，對(duì)于由式(5)代表的峭度與形狀參數(shù)間的函數(shù)K=g(p)及其反函數(shù)p=g-1(K)也無(wú)法給出顯示的表達(dá)式。

由于實(shí)際道路的峭度值取值范圍有限，因此不妨在圖1縱坐標(biāo)取值范圍內(nèi)對(duì)該曲線(xiàn)的反函數(shù)曲線(xiàn)進(jìn)行擬合，得到如下冪函數(shù)形式的p=g-1(K)：

p=2.429K-0.828 8+0.454 1

(6)

圖1 峭度值隨形狀參數(shù)的變化關(guān)系Fig.1 Kurtosis vs.shape factor p

式(6)表明，對(duì)于在區(qū)間[2.83,70]內(nèi)的任意指定的峭度值(該取值范圍已大于實(shí)際道路峭度值的取值范圍)，均可通過(guò)該式快速計(jì)算出所要生成對(duì)稱(chēng)幅度威布爾分布隨機(jī)過(guò)程的形狀參數(shù)p。

當(dāng)已求得形狀參數(shù)p時(shí)，尺度參數(shù)q則可由給定的功率譜密度函數(shù)(曲線(xiàn))確定。

由式(2)可知，Y的二階原點(diǎn)矩為：

(7)

另一方面，當(dāng)給定隨機(jī)過(guò)程Y的功率譜密度PY(f)時(shí)，根據(jù)帕賽瓦爾定理，有：

(8)

當(dāng)給定功率譜密度曲線(xiàn)時(shí)，式(8)中等式右邊積分即代表了曲線(xiàn)與坐標(biāo)軸所包圍的面積。由式(7)、(8)可計(jì)算出尺度參數(shù)：,

(9)

2 指定功率譜密度的對(duì)稱(chēng)幅度威布爾分布隨機(jī)變量產(chǎn)生方法

為了產(chǎn)生對(duì)稱(chēng)幅度威布爾分布隨機(jī)變量，首先考慮產(chǎn)生威布爾分布隨機(jī)變量。

ZMNL法使用兩次非線(xiàn)性變換來(lái)生成指定功率譜密度的威布爾分布隨機(jī)變量，其中第一次非線(xiàn)性變換為通過(guò)目標(biāo)功率譜密度來(lái)(隱式地)計(jì)算高斯隨機(jī)變量的功率譜密度(相關(guān)系數(shù))：

sij=IFFT(PY(ω))=

(10)

式中：2F1(·)為高斯超幾何函數(shù)；sij、ρij分別為威布爾及高斯分布隨機(jī)變量的相關(guān)系數(shù)。

第二次非線(xiàn)性變換由兩組獨(dú)立且各自相關(guān)系數(shù)為式(10)計(jì)算所得的高斯隨機(jī)變量產(chǎn)生威布爾分布隨變量：

(11)

式中：n1、n2～N(0,σ2)為服從高斯分布的隨機(jī)變量，且方差σ2與所要產(chǎn)生的威布爾隨機(jī)變量的尺度參數(shù)的關(guān)系為：

(12)

然而，式(10)所示的非線(xiàn)性變換并不適用于對(duì)稱(chēng)幅度威布爾隨機(jī)變量，且由式(11)所給出的非線(xiàn)性變換總無(wú)法產(chǎn)生負(fù)值變量。因此，ZMNL方法并不能直接用于生成對(duì)稱(chēng)幅度威布爾分布變量。

當(dāng)已由式(6)、(9)分別求得幅度威布爾分布參數(shù)p,q時(shí)，本文保留式(11)所示的第二次非線(xiàn)性變換，并采取如圖2所示的流程來(lái)實(shí)現(xiàn)對(duì)稱(chēng)幅度威布爾分布隨機(jī)變量的生成。

圖2 對(duì)稱(chēng)幅度威布爾分布隨機(jī)變量產(chǎn)生流程Fig.2 Generating process of symmetrical-Weibull-range random variables

圖2過(guò)程中先產(chǎn)生指定功率譜的高斯隨機(jī)變量z(由白噪聲wn3生成指定功率譜的高斯序列z的過(guò)程可由多種現(xiàn)有成熟方法實(shí)現(xiàn)，圖中所示僅為原理性框圖)，將其排序后所得的大小順序?qū)τ墒?11)生成并符號(hào)隨機(jī)化后的對(duì)稱(chēng)幅度威布爾分布隨機(jī)變量x進(jìn)行重新排序。根據(jù)式(9)，源于同一目標(biāo)功率譜的非高斯隨機(jī)變量x與高斯隨機(jī)變量z具有相同的均方根值，加之幅度威布爾分布與高斯分布的相似性，當(dāng)隨機(jī)變量長(zhǎng)度足夠長(zhǎng)時(shí)，在除極值部分的大部分取值范圍內(nèi)，可認(rèn)為由上述排序所得的非高斯變量y是對(duì)z的“仿形”，又由于排序操作不改變幅度分布，因此，重排后的隨機(jī)變量y將同時(shí)具有與x完全一致的幅度分布以及與相似的相關(guān)系數(shù)(即功率譜密度)。

3 試驗(yàn)場(chǎng)道路譜重構(gòu)實(shí)例

圖3所示為某汽車(chē)試驗(yàn)場(chǎng)內(nèi)某可靠性試驗(yàn)道路的功率譜[14]。為驗(yàn)證本文所述方法的有效性，以該功率譜為目標(biāo)進(jìn)行道路總長(zhǎng)為5 km的路形重構(gòu)。

圖3 某試驗(yàn)場(chǎng)可靠性試驗(yàn)道路功率譜密度Fig.3 PSD of durability test road of proving ground

當(dāng)分別指定峭度為4、6、8時(shí)，根據(jù)式(6)及式(9)，所需威布爾分布參數(shù)見(jiàn)表1。

表1 不同峭度下的威布爾分布參數(shù)

為便于比較，分別使用諧波疊加法和本文所述方法對(duì)該試驗(yàn)場(chǎng)路譜進(jìn)行重構(gòu)。重構(gòu)路形譜及其功率譜密度分別如圖4、5所示。

圖4 不同峭度下的重構(gòu)路形Fig.4 Reconstructed road profiles with different kurtosis

圖5 目標(biāo)功率譜與重構(gòu)功率譜Fig.5 Target PSD vs. reconstructed PSDs

圖6 重構(gòu)路形局部放大圖Fig.6 Local details of reconstructed road profile

將圖4所示高斯路形及K=4時(shí)的非高斯路形繪制在同一圖中并局部放大，如圖6所示。由該局部路形可知由本文方法生成的路形與高斯路形在大體相似的同時(shí)，具有高處更高、低處更低的非高斯特性。

不同峭度下的重構(gòu)非高斯路形與高斯路形譜的幅值概率密度分布如圖7所示。由圖7可知，當(dāng)峭度值較為接近3時(shí)(此時(shí)形狀參數(shù)p大于1)，對(duì)稱(chēng)威布爾分布將在接近0幅值附近出現(xiàn)對(duì)稱(chēng)的雙峰，當(dāng)峭度增大時(shí)，由于威布爾分布退化為指數(shù)分布，則其呈雙指數(shù)分布。因此，本文所述方法對(duì)于描述那些路況惡劣、峭度值較高的道路效果更為理想。

上述重構(gòu)結(jié)果表明，本文方法生成的非高斯相路形比傳統(tǒng)高斯路形，可以在保持功率譜的同時(shí)，提供更強(qiáng)烈的路面激勵(lì)，且明確的幅值概率密度分布可為后續(xù)試驗(yàn)時(shí)進(jìn)行載荷處理提供依據(jù)，驗(yàn)證了本文所述方法的正確性、實(shí)用性和有效性。同時(shí)，該結(jié)果也直觀地表明了以往單純依靠功率譜來(lái)描述道路激勵(lì)特征的不足之處，說(shuō)明了進(jìn)行非高斯路譜重構(gòu)的必要性。

(a) K=4

(b) K=6

(c) K=8圖7 不同峭度下的重構(gòu)路形的幅值概率密度Fig.7 PDFs of road profiles with different kurtosis

4 結(jié) 論

本文由威布爾分布的高階矩特性推導(dǎo)了具有指定峭度的對(duì)稱(chēng)幅度威布爾分布隨機(jī)過(guò)程分布參數(shù)的計(jì)算公式；以ZMNL方法為基礎(chǔ)，提出了由高斯白噪聲產(chǎn)生指定峭度值與功率譜密度對(duì)稱(chēng)幅度威布爾分布隨機(jī)變量的方法。通過(guò)將本文方法應(yīng)用于某試驗(yàn)場(chǎng)可靠性試驗(yàn)道路譜重構(gòu)，實(shí)現(xiàn)了產(chǎn)生具有指定峭度及功率譜的非高斯路形。與現(xiàn)有方法產(chǎn)生的高斯路形進(jìn)行比較，證明了該非高斯路形更適用于具有較高嚴(yán)酷程度的各種試驗(yàn)道路的重構(gòu)，驗(yàn)證了本文所述方法的實(shí)用性和有效性。

由于本文所述方法對(duì)于功率譜及重構(gòu)隨機(jī)變量的單位不做限制，因此該方法亦可廣泛應(yīng)用于對(duì)各類(lèi)隨機(jī)振動(dòng)試驗(yàn)所需的載荷譜重構(gòu)工作。

[ 1 ] DODDS C J, ROBSON J D. The description of road surface roughness [J]. Journal of Sound and Vibration, 1973, 1(2): 175-183.

[ 2 ] 檀潤(rùn)華, 陳鷹, 路甬祥. 路面對(duì)汽車(chē)激勵(lì)的時(shí)域模型建立及計(jì)算機(jī)仿真[J]. 中國(guó)公路學(xué)報(bào), 1998, 11(3): 96-101.

TAN Runhua, CHEN Ying, LU Yongxiang. The mathematical models in time domain for the road disterbances and the simulation [J]. China Journal of Highway and Transport, 1998, 11(3): 96-101.

[ 3 ] 劉宏偉, 陳綿書(shū), 王勛龍. 路面譜室內(nèi)再現(xiàn)方法的研究[J]. 實(shí)用測(cè)試技術(shù), 1999, 25(5): 17-18.

LIU Hongwei, CHEN Mianshu, WANG Xunlong. Study on road profile indoor reconstruction [J]. Practical Measurement Technology, 1999, 25(5): 17-18.

[ 4 ] 張永林, 鐘毅芳. 車(chē)輛路面不平度輸入的隨機(jī)激勵(lì)時(shí)域模型 [J]. 農(nóng)業(yè)機(jī)械學(xué)報(bào), 2004, 35(2): 9-12.

ZHANG Yonglin, ZHONG Yifang. Time domain model of road undulation excitation to vehicles [J]. Transactions of The Chinese Society of Agricultural Machinery, 2004, 35(2): 9-12.

[ 5 ] 張永林. 用諧波疊加法重構(gòu)隨機(jī)道路不平順高程的時(shí)域模型 [J]. 農(nóng)業(yè)工程學(xué)報(bào), 2003, 19(6): 32-34.

ZHANG Yonglin. Time domain model of road irregularities simulated using the harmony superposition method [J]. Transactions of The Chinese Society of Agricultural Engineering, 2003, 19(6): 32-34.

[ 6 ] BOGSJ? K. Road profile statistics relevant for vehicle fatigue [D]. Lund,Sweden: Lund University, 2007.

[ 7 ] 叢楠, 尚建忠, 任焱晞. 基于對(duì)稱(chēng)α-穩(wěn)定分布的可靠性試驗(yàn)場(chǎng)道路時(shí)域激勵(lì)建模方法 [J]. 汽車(chē)工程, 2012, 24(12): 32-34.

CONG Nan, SHANG Jianzhong, REN Yanxi. Modeling methodology of time domain reliability test ground road disturbances based on symmetrical α-stable distribution [J]. Automotive Engineering, 2012, 24(12): 32-34.

[ 8 ] 劉凡, 艾加秋. 用ZMNL方法實(shí)現(xiàn)地面雜波的建模與仿真 [J]. 中國(guó)科學(xué)院研究生院學(xué)報(bào), 2010, 27(2): 275-279.

LIU Fan, AI Jiaqiu. Modeling and simulation of ground clutter using ZMNL algorithm [J].Journal of Graduate School of the Chinese Academy of Science, 2010, 27(2): 275-279.

[ 9 ] 李青華, 孔令講, 楊曉波. 基于SIRP法的相關(guān)韋布爾分布雷達(dá)雜波仿真[J]. 雷達(dá)科學(xué)與技術(shù), 2011, 9(3): 253-258.

LI Qinghua, KONG Lingjiang, YANG Xiaobo. Simulation of correlated weibull distribution radar clutter based on sirp method [J].Radar Science and Technology, 2011, 9(3):253-258.

[10] SMALLWOOD D O. Generation of time histories with a specified auto spectral density and probability density function[J]. Shock and Vibration, 1997, 4(5): 361-377.

[11] 蔣瑜, 陳循, 陶俊勇, 等. 指定功率譜密度、偏斜度和峭度值下的非高斯隨機(jī)過(guò)程數(shù)字模擬 [J]. 系統(tǒng)仿真學(xué)報(bào), 2006, 18(5): 1127-1130.

JIANG Yu, CHEN Xun, TAO Junyong, et al. Numerically simulationn non-gaussian random preocesses with specified PSD, skewness and kurtosis [J]. Journal of System Simulation, 2006, 18(5): 1127-1130.

[12] 陳金明. 雷達(dá)雜波建模與仿真[D]. 西安：西安電子科技大學(xué), 2010.

[13] 姜培華. 幾種概率分布高階原點(diǎn)矩的計(jì)算 [J]. 重慶工商大學(xué)學(xué)報(bào), 2014,31(9): 1-5.

JIANG Peihua. Calculation of higher-order origin moment of several probability distributions [J]. Journal of Chongqing Technology and Business University, 2014,31(9): 1-5.

[14] 虞明, 趙濟(jì)海, 鄔惠樂(lè), 等. 定遠(yuǎn)試驗(yàn)區(qū)汽車(chē)試驗(yàn)路面譜的研究 [J]. 汽車(chē)工程, 1990, 49(2): 22-30.

YU Ming, ZHAO Jihai, WU Huile, et al. Road spectral analysis of dingyuan vehicle testing ground [J]. Automotive Engineering, 1990, 49(2): 22-30.