路面功率譜密度換算及不平度建模理論研究*

杜 峰， 葛曉成， 陳 翔， 丁金全

(1.華南理工大學機械與汽車工程學院 廣州，510641) (2.重慶大學機械傳動國家重點實驗室 重慶，400030) (3.吉林大學汽車仿真與控制國家重點實驗室 長春，130022)

?

路面功率譜密度換算及不平度建模理論研究*

杜 峰1， 葛曉成2， 陳 翔3， 丁金全3

(1.華南理工大學機械與汽車工程學院 廣州，510641) (2.重慶大學機械傳動國家重點實驗室 重慶，400030) (3.吉林大學汽車仿真與控制國家重點實驗室 長春，130022)

為了明晰路面不平度空間域統計量的計算，以及幾種重要功率譜密度(power spectral density,簡稱PSD)之間的關系，以帕塞瓦爾定理和維納-辛欽定理為依據，在推導空間域自相關函數和PSD計算公式的基礎上，導出了不平度空間域位移、速度與加速度PSD以及空間頻率與角頻率PSD之間的換算關系。另外，為了完善傅里葉逆變換法路面建模中PSD離散化的理論基礎，以傅里葉級數與變換、離散傅里葉變換和頻域卷積定理為依據，從離散化的原因、目的和結果驗證出發論證了PSD離散化的正確性。就模擬路面驗證問題，指出直接法譜估計的不合理之處，論證了平均周期圖法譜估計時，空間與時間采樣頻率分別對應著空間域和時間域PSD輸出。結果表明，上述換算關系和論證是正確的，可應用于路面不平度建模和汽車平順性分析。

公路； 汽車平順性； 功率譜密度； 路面不平度建模

引 言

空間域隨機信號分析對于路面不平度建模和汽車平順性分析非常重要[1-3]，大多文獻通過車速將其轉換到時間域進行分析，影響了對路面空間域統計量的計算和理解。

空間域路面不平度位移、速度和加速度功率譜密度之間的換算關系，是汽車平順性分析的重要基礎。文獻[4]已給出確切表達式，但沒有解釋原由。文獻[5]也沒給出推導，現有文獻大多直接應用其結論[6]。盧士富[7]研究了這個問題，現在加以推導。

路面不平度空間頻率功率譜密度Gq(n)和空間角頻率功率譜密度Gq(Ω)之間的數量關系，在運用維納-辛欽定理，經由功率譜密度計算隨機信號的標準差時非常重要，模擬路面標準差是驗證生成路面準確性的基本依據之一，但對此關系的分析論證和檢驗很少，散布在文獻中的應用顯示對二者關系理解不夠全面，需要做深入研究。

針對上述現象，筆者將以帕塞瓦爾定理和維納-辛欽定理為理論依據，對路面空間域自相關函數和功率譜密度的計算做出分析，并推導出路面不平度位移、速度和加速度功率譜密度之間換算關系，根據傅里葉變換對的兩種形式和文獻[4-5]，論證了空間頻率和角頻率功率譜密度之間的兩種對應關系。

在汽車平順性仿真中，基于功率譜密度離散化的傅里葉逆變換法是路面空間建模的常用方法之一，其建模思想比較成熟，但尚有一個關鍵步驟缺乏嚴格理論分析[8]——功率譜密度離散化。現有文獻通常都略過這些問題[9]。

另外，在應用傅里葉逆變換法做路面模擬時，細節處理不當會導致錯誤。奇偶采樣點的對稱設置有區別，功率譜密度驗證不能采用直接法計算，平均周期圖法譜估計采用空間和時間采樣頻率分別對應著不同的功率譜密度輸出，這類關鍵細節在之前文獻很少涉及，文中將逐一分析。

筆者在深入理解傅里葉變換不同分類(周期連續信號傅里葉級數和傅里葉變換、能量信號傅里葉變換、離散時間傅里葉變換、離散傅里葉變換和周期離散信號傅里葉級數)本質特性的基礎上，綜合運用傅里葉變換思想和頻域卷積定理，從理論上嚴格分析了功率譜密度離散化的問題，并對工程應用中的幾個關鍵細節給予分析，完善了基于功率譜密度離散化的傅里葉逆變換法路面建模的理論基礎，為準確模擬路面不平度，進行汽車平順性分析提供了保障。

1 空間域路面不平度計算

頻率可以指空間、時間和其它虛擬域，是單位域長度所包含波長的個數。路面建模在空間域中研究路面形狀，故采用空間頻率n，指單位空間(1 m)中所包含的路面波形的波長個數。空間角頻率Ω=2πn。空間頻率和角頻率分別對應時間域中的f和2πf，空間域信號變量l對應時間域變量t。

1.1 位移不平度自相關函數和功率譜密度

路面位移不平度q=q(l)，其中自變量是長度，因變量是高度，均為長度單位量綱。q(l)為時域無限信號，工程中通常截取一段長度L進行研究。L的取值，要保障空間頻率分辨率，對于路面不平度信號的空間頻率在[0.011 2.83]m-1之間，故最小辨識頻率dn=1/L≤0.011，即L≥91 m。

空間長度L的路面不平度信號為

(1)

其自相關函數為

(2)

其中：qL(l)為空間域有限信號，其空間域頻譜為連續譜Q(n,L)。

根據帕塞瓦爾定理，路面不平度信號的平均功率為

(3)

空間域單邊功率譜密度Gq(n)為

(4)

由維納-辛欽定理可換算出路面位移不平度時間域功率譜密度Gq(f)。

(5)

其中：l=ut(位移等于車速與時間乘積)，f=nu(時間頻率等于空間頻率與車速乘積)。

1.2 位移和速度及加速度功率譜密度關系

由自相關定義式(2)可得

(6)

自相關函數對長度變量l的一、二階導數[7]分別為

(7)

(8)

根據維納-辛欽定理，改變求導和積分次序得

(9)

(10)

由式(8～10)可得

(11)

同理將自相關函數求三、四階導數，可得

(12)

1.3 空間頻率和空間角頻率功率譜密度關系

首先將傅里葉變換的常用形式表述如下

(13)

GB/T 7031—2005和ISO 8608—1995中的一個重要關系式(14)可以作為空間頻率功率譜密度Gq(n)與角頻率功率譜密度Gq(Ω)關系論證的依據。

Gq(n0)=16Gq(Ω0)

(14)

統計分析表明，路面隨機高程信號滿足零均值高斯分布，即q～N(0,σ2)，其均方值與零位置自相關函數值相等。

(15)

根據維納-辛欽定理，由式(13)可知零位置自相關函數值為

(16)

由式(15，16)可知

Sq(Ω)=Sq(n)

(17)

根據式(14)做一驗證，雖然Ω=2πn，但不包括Ω0和n0點，其值[4-5]分別為1 rad/m和0.1 m-1。

(18)

(19)

由式(17～19)可得，Gq(n0)=2.533Gq(Ω0)，與式(14)不同。注意到傅里葉變換對的另一種表達，如式(20)所示，則自相關函數如式(21)所示。

(20)

(21)

由式(15)，(21)得到Gq(Ω)dΩ=Gq(n)dn，故

Sq(n)=2πSq(Ω)

(22)

將式(18)，(19)代入式(22)，得Gq(n0)=16Gq(Ω0)，與國標一致，說明在國標和ISO文件中，由采樣路面信號自相關函數求路面位移功率譜密度時，采用的傅里葉變換為式(20)，而不是通常采用的式(13)，同時也驗證了關系式(22)的正確性。

綜合上述分析，空間頻率和角頻率功率譜密度的關系，隨著維納-辛欽定理使用時所采用的傅里葉變換對的形式而不同，即關系式(17)對應式(13)，關系式(22)對應式(20)。

關系式(17)對應式(13)的正確性檢驗將在1.4.1節，對應用這個數量關系計算的路面不平度位移標準差作校驗時進行；關系式(22)對應式(20)正確的原因在于，它能導出Gq(n0)=16Gq(Ω0)，與國標一致。

1.4 空間域路面不平度位移標準差計算

(23)

統計顯示路面空間頻率分布于[0.011 2.83] m-1區間,在常用車速下(10～30 m/s)下，激勵的時間頻率(車速與空間頻率乘積)分別為[0.11 28.3]和[0.33 84.9] Hz，兩區間交集為[0.33 28.3] Hz，涵蓋了車身(1～2 Hz)和車輪(10～20 Hz)的固有頻率，適合做汽車平順性分析的輸入。根據國家標準，A級路面不平度系數幾何均值為16×10-6m3，則其不平度路面標準差為

(24)

驗證依據GB/T 7031-2005，國標附錄C中表C.3列出不同等級道路位移功率譜密度的幾何均值和限值，對于A級道路，對應于空間中心頻率為0.015 6，0.031 2，0.062 5，0.125，0.25，0.5，1和2 m-1的倍頻段的單邊位移功率譜密度幾何均值分別為655×10-6，164×10-6，41×10-6，10.2×10-6，2.56×10-6，0.64×10-6，0.16×10-6和0.04×10-6m3。

中心頻率點0.015 6，2 m-1對應的倍頻段分別為[0.011 0 0.022 1]，[1.414 2 2.828 4] m-1，可見，自0.015 6到2 m-1的倍頻帶剛好涵蓋了路面信號的頻率范圍[0.011 2.83] m-1，對其中各倍頻段位移功率譜密度幾何均值求和開方可求得A級路面平均標準差，約等于3.78 mm。

可見式(24)的計算結果與標準值近似相等，從而驗證了空間頻率和角頻率功率譜密度的關系(即式(17)與式(13)對應)的正確性。

另外，如果在時間域計算標準差，如式(25)，(26)所示，結果與式(24)相同，說明路面均方根與車速和計算所在域無關。

(25)

(26)

2 傅里葉逆變換法路面模擬

要從頻域數據合成路面信號，必須知道幅值譜和相位譜。統計分析發現路面相位譜服從[0 2π]的均勻分布，可由函數隨機生成。基于功率譜密度離散化的路面建模方法，就是要由已知的標準等級路面功率譜密度計算出信號的幅值譜，進而由傅里葉逆變換法得到路面不平度信號。

傅里葉變換有多種類型：周期連續信號傅里葉級數(FS)、周期連續信號傅里葉變換、非周期連續信號傅里葉變換(FT)、離散時間非周期信號傅里葉變換、離散傅里葉變換(DFT)和周期離散信號傅里葉級數(DTFS)。只有DFT和DTFS在時、頻域均離散，并且二者定義完全相同，DTFS是DFT的周期延拓，計算機只做DFT。2.2節將證明DFT變換結果可看作信號幅值譜，用于路面信號合成。

2.1 功率譜密度離散化的原因及目的

路面不平度可看作各態歷經隨機信號，取長度為L的一段路面qL(l)，根據式(4)可知，路面信號的單邊功率譜密度估計為

(27)

qL(l)為空間域有限的非周期連續信號，其傅里葉變換是連續值，分布在整個空間頻域，理論上其幅值譜(相當于周期信號傅里葉級數的傅里葉系數)趨于零，由式(4)和(27)知Q(n,L)是信號的幅值譜密度，而信號合成需要的是幅值譜Q(k)，故不能直接用來合成信號。

雖然理論上其幅值譜趨于零，但有效路面數據只是[0L]這一段，將其周期延拓成為周期為L的周期信號，滿足狄利克雷條件，傅里葉級數存在，可以由傅里葉級數(即幅值譜)合成一個周期[0L]的路面不平度信號，如式(28)所示，這里Ω0=2π/L為空間角頻率(基頻)。

(28)

為了適應計算機處理，將周期延拓的路面不平度信號進行等周期采樣，則信號成為周期離散時間信號，于是其時、頻域均離散，從而可以由DFT的定義式(29)求出路面不平度信號的幅值譜Q(k)。

(29)

將Q(n,L)積分形式表示成N項級數和，與離散傅里葉變換Q(k)產生換算關系。

大多數商圈只注重交易數據而忽視購買群體的過程數據，而這部分數據正是挖掘商機的入口。圖1為商圈客流分析功能示意,其中客流分析項是對商圈日常的經營數據和客流數據進行匯總得到的，這些數據的處理和分析對于商圈運營者而言可進行更為細致的分析。客流數據分析系統需結合交易數據和過程數據獲取有價值的數據。研究結果證明，只有0.5%～3.0%的人在網上點擊廣告之后會真正購買相應的產品,而消費者走進實體商店購買的概率則大很多。對于商圈管理者來說，該部分數據的價值分析示意見圖2。

(30)

其中：Δn=1/(NΔl)，為空間頻率分辨率。

從以上分析可看到，功率譜密度離散化本質是將截斷路面信號的幅值譜密度離散化，得到延拓周期信號的傅里葉級數，以此來合成路面信號。

2.2 離散化結果驗證

功率譜密度離散化的目的是得到信號的幅值譜，那么Q(k)是不是路面不平度離散周期信號的幅值譜？下面從周期連續信號傅里葉級數(式(28))出發來給以驗證。

將長度和周期均為L的路面不平度連續信號的傅里葉級數近似表達成M項級數和的形式，可得到連續周期信號傅里葉級數與DFT的關系

(31)

可以看出離散傅里葉變換Q(k)與幅值級數Ck只存在一個系數差別，這個系數將在逆變換時得以消除，故Q(k)可以代表信號幅值譜。

亦可如此分析,假設信號qL(l)長度和周期L均為無窮大，根據幅值譜和幅值譜密度的關系得到

(32)

對路面qL(l)在一個周期內等間隔采樣N點，根據頻域卷積定理，信號頻譜被延拓成N個，路面客觀存在，即傅里葉級數不變，于是得到

(33)

(34)

比較式(30)和式(34)發現，Q(k)即為信號幅值譜。采樣點數N的取值可以參考文獻[9]。

2.3 頻譜點對稱設置

根據上述分析，離散功率譜密度估計式為

(35)

信號幅值譜為

(36)

其中：dn=1/L，為空間頻率分辨率，即采樣路面總長度的倒數。

相位角在[0 2π]上均勻分布，在Matlab中可以由unifrnd函數生成。然后根據式(36)通過傅里葉逆變換生成路面空間域隨機信號。

離散傅里葉變換Q(k)的第一個點(直流分量)不存在對稱。對于模擬路面數據N，如果是偶數，則關于N/2+1點對稱；若為奇數，則關于(N+2)/2點對稱，根據式(36)生成幅值譜Q(k) (k=1,…,N/2+1或(N+1)/2)。則另一半幅值頻譜數據為

Q*(N/2)，Q*(N/2-1)，…，Q*(2) (N為偶數)；Q*((N+1)/2),…,Q*(2) (N為奇數)。

Q*(k)表示Q(k)的共軛復數，將兩段頻譜數據組合后進行離散傅里葉逆變換即可得到路面信號。

3 功率譜密度驗證

正確的模擬路面除了要滿足均值和標準差的條件外(即均值接近于零；標準差與標準值接近并不隨車速、采樣頻率和仿真時間改變)，功率譜密度曲線的趨勢必須與國標和ISO文件要求一致[10]。

3.1 功率譜密度驗證方法

計算模擬路面的功率譜密度時，不少文獻采用了直接法[9]，這是不準確的。因為直接法采用的計算式為式(35)，本質上是上述計算的逆過程，這樣得到的功率譜密度數值與標準值完全相同，如圖1所示，模擬路面的譜曲線與ISO標準路面譜曲線完全重合，驗證沒有意義，不能真實反映程序模擬路面的頻域特性。

圖1 直接法功率譜密度計算Fig.1 PSD calculated by direct computation method

經典譜估計的平均周期圖Welch法，對路面信號分段，段間適當重疊，分段計算功率譜密度再平均，通過窗函數減小信號分段造成的能量向旁瓣泄漏，窗函數造成的信號能量衰減，已在函數內部進行了修正，Welch法能更真實地反映模擬路面信號功率在頻率軸的分布特征。

3.2 兩種采樣頻率的應用

平均周期圖法譜估計函數pwelch調用格式為：[Pxx,F]=pwelch(signal,window,overlap,NFFT,fs,'oneside')，有兩種采樣頻率：空間域采樣頻率ns和時間域采樣頻率fs，令車速為u，采樣時間間隔為Δt，則時間和空間頻率分別為：fs=1/Δt，ns=1/(uΔt)。

窗函數選擇，國標推薦用漢寧窗，這其實是在頻率分辨率和能量泄漏之間取的一個較理想折中方案，這里靈活采用切比雪夫窗，旁瓣衰減100 dB，能量泄漏減小，但頻率分辨力下降，譜曲線趨于平滑。

窗函數通常應用于濾波器設計和譜估計，它們對窗函數的要求有差異：譜估計防止能量泄漏，希望窗函數頻譜旁瓣較小；濾波器設計通常要求過渡帶較小，要求窗函數主瓣陡峭，但有關指標是矛盾的，需要折中，詳細的選擇方法可參考有關教材。

采用空間域采樣頻率ns將得到空間域路面不平度位移功率譜密度Gq(n)，同時輸出的頻率點是空間頻率；如果采用時間域采樣頻率fs，則得到時間域的頻率點和相應的時域功率譜密度Gq(f)。

3.2.1 空間域路面不平度位移功率譜密度

應用Welch法對模擬路面信號Xm進行單邊功率譜密度估計，根據上述分析，編制Matlab程序如下，運行結果如圖2所示。

NFFT=2^nextpow2(round(length(Xm)/4));

window=chebwin(NFFT); %切比雪夫窗函數

overlap=round(length(Xm)/8); %段間重疊

ns=ceil(N/L); %空間采樣頻率

[Pxx,F]=pwelch(Xm,window,overlap,NFFT,ns,'oneside'); %采用空間采樣頻率做功率譜密度估計

圖2 采用空間采樣頻率的Welch法輸出Fig.2 Output of Welch with spatial sampling frequency

圖2中虛線為B級路面位移功率譜密度函數曲線，數學函數為Gq(n)=64×10-6×(n/0.1)-2，在雙對數坐標中為一直線；程序中ns為空間采樣頻率，圖中實線為模擬路面功率譜密度曲線，在標準曲線上下波動，趨勢一致，說明采用空間采樣頻率得到的是空間域功率譜密度。

3.2.2 時域路面不平度位移功率譜密度

采用時間采樣頻率，將程序參數ns改為時域采樣頻率fs，運行程序得到圖3。

圖3 采用時間采樣頻率的Welch法輸出Fig.3 Output of Welch with temporal sampling frequency

車速u為20 m/s，時間頻率f=20n，即輸出的頻率序列F為時域點，故圖2和圖3的橫坐標不同。圖3的虛線數學函數為Gq(f)=Gq(n)/u=3.2×10-6×(n/0.1)-2，為車速20 m/s下B級標準路面的時域功率譜密度，圖3顯示模擬路面的功率譜密度與標準時域功率譜密度一致，故可知pwelch函數采用時間采樣頻率時輸出的功率譜密度和頻率點均為時域。

如果要換算成空間域的功率譜密度，根據空間與時間域頻率和功率譜密度的關系，需將輸出的頻率向量F除以車速u，相應的功率譜密度Pxx乘以車速。

4 結束語

以帕塞瓦爾定理和維納-辛欽定理為基礎，分析了路面空間域自相關函數和功率譜密度的計算方法，推導出路面不平度空間域位移、速度和加速度功率譜密度之間，以及空間頻率與空間角頻率功率譜密度之間換算關系。應用傅里葉級數和變換、離散傅里葉變換、離散時間傅里葉級數理論和頻域卷積定理，詳細分析論證了傅里葉逆變換路面建模方法的功率譜密度離散化問題。指出了直接法功率譜密度驗證的不合理之處，論證了平均周期圖法譜估計時應用空間采樣頻率將輸出路面空間域PSD，而時間采樣頻率對應著路面時間域PSD輸出。理論分析以及與國標和ISO標準的對比驗證顯示，文中的研究是正確的，可用于路面不平度建模和汽車平順性分析。

[1] Bogsj? K, Podgórski K, Rychlik I. Models for road surface roughness [J]. Vehicle System Dynamics, 2012, 50(5): 725-747.

[2] Mucka P. Influence of road profile obstacles on road unevenness indicators [J]. Road Materials and Pavement Design, 2013, 14(3): 689-702.

[3] 黃雪濤，顧亮，呂唯唯. 輕型載貨汽車減振技術分析及優化分析[J]. 振動、測試與診斷，2013, 33(2): 315-318.

Huang Xuetao, Gu Liang, Lv Weiwei. Vibration reduction technology analysis and optimal design on lighttruck [J]. Journal of Vibration, Measurement & Diagnosis, 2013, 33(2): 315-318. (in Chinese)

[4] International Organization for Standardization. ISO 8608—1995(E), Mechanical vibration-road surface profiles-reporting of measured data [S]. Geneva: Int. Standard Organization, 1995.

[5] 全國機械振動與沖擊標準化技術委員會. GB/T 7031—2005 機械振動——道路路面譜測量數據報告[S].北京:中國標準出版社,2006.

[6] Agostinacchio M, Cinampa D, Olita S. The vibrations induced by surface irregularities in road pavements-a Matlab?approach [J]. European Transport Research Review, 2014, 6(3): 267-275.

[7] 盧士富. 關于平穩過程自譜密度運算中的兩個問題[J]. 西安公路學院學報，1991,11(2): 54-61.

Lu Shifu. On two problems in calculating self-spectrum density of smooth process [J]. Journal of Xi′an University of Highway, 1991, 11(2): 54-61. (in Chinese)

[8] Zhang Yonglin, Zhang Jiafan. Numerical simulation of stochastic road process using white noise filtration [J]. Mechanical Systems and Signal Processing, 2006, 20(2): 363-372.

[9] 劉獻棟，鄧志黨，高峰. 基于逆變換的路面不平度仿真研究[J]. 中國公路學報，2005,18(1): 122-126.

Liu Xiandong, Deng Zhidang, Gao Feng. Study of simulation of road roughness based on inverse transform [J]. China Journal of Highway and Transport, 2005, 18(1): 122-126. (in Chinese)

[10]Durst P J, Mason G L, McKinley B, et al. Predicting RMS surface roughness using fractal dimension and PSD parameters [J]. Journal of Terramechanics, 2011, 48(2): 105-111.

10.16450/j.cnki.issn.1004-6801.2015.05.030

2014-09-24；

2014-12-01

TH113.1; U461.5+6; TP391.9

杜峰，男，1980年5月生，博士生。主要研究方向為汽車平順性仿真及客觀評價。曾發表《PSD和PWELCH函數的分析改進及應用》(《中國測試》2010年第36卷第1期)等論文。 E-mail：dufeng123dufeng@126.com