一种指定功率谱密度与峭度值的对称威布尔分布道路谱重构方法

丛 楠， 陈俊达， 任焱晞， 王彬星， 李青霞

(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.