一種指定功率譜密度與峭度值的對稱威布爾分布道路譜重構方法

叢 楠， 陳俊達， 任焱晞， 王彬星， 李青霞

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

功率譜密度作為路面不平度在(空間)頻域上能量分布的指征，自20世紀70年代起，已成為描述道路不平度的最主要形式[1]。當前，無論是普通公路還是各種專用試驗道路，功率譜密度均是道路不平度的最主要數值特征。在進行道路模擬時，定義在(空間)頻域上的功率譜密度須“重構”為隨機時域或空間域歷程，方能被輸入振動臺架開展試驗。目前成熟的以諧波疊加法和線性濾波為代表的路譜重構方法均僅能產生高斯分布的隨機歷程[2-5]。實測數據表明，高斯路形譜僅適用于路況良好的等級道路，對于路況較差或專用可靠性試驗道路，由于其往往包含有超過“標準”預期的不規則性，使得其幅值分布具有“厚尾”特征而呈現出明顯的非高斯分布特性[6-7]。

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

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

1 指定峭度的對稱幅度威布爾分布隨機過程

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

對于滿足雙參數威布爾分布隨機過程X～Wweibull(p,q)，其具有如式(1)所示的概率密度函數：

(1)

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

根據威布爾分布性質，其k階原點矩為[13]：

(2)

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

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

(3)

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

(4)

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

(5)

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

由于實際道路的峭度值取值范圍有限，因此不妨在圖1縱坐標取值范圍內對該曲線的反函數曲線進行擬合，得到如下冪函數形式的p=g-1(K)：

p=2.429K-0.828 8+0.454 1

(6)

圖1 峭度值隨形狀參數的變化關系Fig.1 Kurtosis vs.shape factor p

式(6)表明，對于在區間[2.83,70]內的任意指定的峭度值(該取值范圍已大于實際道路峭度值的取值范圍)，均可通過該式快速計算出所要生成對稱幅度威布爾分布隨機過程的形狀參數p。

當已求得形狀參數p時，尺度參數q則可由給定的功率譜密度函數(曲線)確定。

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

(7)

另一方面，當給定隨機過程Y的功率譜密度PY(f)時，根據帕賽瓦爾定理，有：

(8)

當給定功率譜密度曲線時，式(8)中等式右邊積分即代表了曲線與坐標軸所包圍的面積。由式(7)、(8)可計算出尺度參數：,

(9)

2 指定功率譜密度的對稱幅度威布爾分布隨機變量產生方法

為了產生對稱幅度威布爾分布隨機變量，首先考慮產生威布爾分布隨機變量。

ZMNL法使用兩次非線性變換來生成指定功率譜密度的威布爾分布隨機變量，其中第一次非線性變換為通過目標功率譜密度來(隱式地)計算高斯隨機變量的功率譜密度(相關系數)：

sij=IFFT(PY(ω))=

(10)

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

第二次非線性變換由兩組獨立且各自相關系數為式(10)計算所得的高斯隨機變量產生威布爾分布隨變量：

(11)

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

(12)

然而，式(10)所示的非線性變換并不適用于對稱幅度威布爾隨機變量，且由式(11)所給出的非線性變換總無法產生負值變量。因此，ZMNL方法并不能直接用于生成對稱幅度威布爾分布變量。

當已由式(6)、(9)分別求得幅度威布爾分布參數p,q時，本文保留式(11)所示的第二次非線性變換，并采取如圖2所示的流程來實現對稱幅度威布爾分布隨機變量的生成。

圖2 對稱幅度威布爾分布隨機變量產生流程Fig.2 Generating process of symmetrical-Weibull-range random variables

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

3 試驗場道路譜重構實例

圖3所示為某汽車試驗場內某可靠性試驗道路的功率譜[14]。為驗證本文所述方法的有效性，以該功率譜為目標進行道路總長為5 km的路形重構。

圖3 某試驗場可靠性試驗道路功率譜密度Fig.3 PSD of durability test road of proving ground

當分別指定峭度為4、6、8時，根據式(6)及式(9)，所需威布爾分布參數見表1。

表1 不同峭度下的威布爾分布參數

為便于比較，分別使用諧波疊加法和本文所述方法對該試驗場路譜進行重構。重構路形譜及其功率譜密度分別如圖4、5所示。

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

圖5 目標功率譜與重構功率譜Fig.5 Target PSD vs. reconstructed PSDs

圖6 重構路形局部放大圖Fig.6 Local details of reconstructed road profile

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

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

上述重構結果表明，本文方法生成的非高斯相路形比傳統高斯路形，可以在保持功率譜的同時，提供更強烈的路面激勵，且明確的幅值概率密度分布可為后續試驗時進行載荷處理提供依據，驗證了本文所述方法的正確性、實用性和有效性。同時，該結果也直觀地表明了以往單純依靠功率譜來描述道路激勵特征的不足之處，說明了進行非高斯路譜重構的必要性。

(a) K=4

(b) K=6

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

4 結 論

本文由威布爾分布的高階矩特性推導了具有指定峭度的對稱幅度威布爾分布隨機過程分布參數的計算公式；以ZMNL方法為基礎，提出了由高斯白噪聲產生指定峭度值與功率譜密度對稱幅度威布爾分布隨機變量的方法。通過將本文方法應用于某試驗場可靠性試驗道路譜重構，實現了產生具有指定峭度及功率譜的非高斯路形。與現有方法產生的高斯路形進行比較，證明了該非高斯路形更適用于具有較高嚴酷程度的各種試驗道路的重構，驗證了本文所述方法的實用性和有效性。

由于本文所述方法對于功率譜及重構隨機變量的單位不做限制，因此該方法亦可廣泛應用于對各類隨機振動試驗所需的載荷譜重構工作。

[ 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 ] 檀潤華, 陳鷹, 路甬祥. 路面對汽車激勵的時域模型建立及計算機仿真[J]. 中國公路學報, 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 ] 劉宏偉, 陳綿書, 王勛龍. 路面譜室內再現方法的研究[J]. 實用測試技術, 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 ] 張永林, 鐘毅芳. 車輛路面不平度輸入的隨機激勵時域模型 [J]. 農業機械學報, 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 ] 張永林. 用諧波疊加法重構隨機道路不平順高程的時域模型 [J]. 農業工程學報, 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 ] 叢楠, 尚建忠, 任焱晞. 基于對稱α-穩定分布的可靠性試驗場道路時域激勵建模方法 [J]. 汽車工程, 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方法實現地面雜波的建模與仿真 [J]. 中國科學院研究生院學報, 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法的相關韋布爾分布雷達雜波仿真[J]. 雷達科學與技術, 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]. 系統仿真學報, 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]. 西安：西安電子科技大學, 2010.

[13] 姜培華. 幾種概率分布高階原點矩的計算 [J]. 重慶工商大學學報, 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]. 汽車工程, 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.