利用聲信號基于HHT方法識別螺旋槳模態參數

劉正浩， 丁 舉， 毛獻群

(1. 中國船舶及海洋工程設計研究院， 上海 200011； 2. 噴水推進技術重點實驗室， 上海 200011)

隨著對螺旋槳靜音性能的要求愈加嚴格，在設計中對螺旋槳的動力響應特性更加重視[1-2]。螺旋槳的模態參數是進行動力響應分析的前提和基礎[3]，在這方面，黃政等[4-7]開展了螺旋槳模態的數值與試驗研究，這些研究中，分別是通過加速度傳感器和電阻應變片對槳葉動態響應進行測量，得到固有頻率和模態阻尼比等參數。吳武輝等還進一步研究了傳感器質量效應對于頻率測量結果的影響，說明隨著傳感器數量增加，試驗得到的固有頻率有明顯的降低。實際上，螺旋槳模態試驗常常針對縮尺模型展開，對于縮尺槳模，單槳葉和傳感器的質量對比并未達到可忽略傳感器質量效應的程度，并且，高階模態由于模態有效質量更小，對于傳感器的質量效應更加敏感，試驗結果可能帶有更大誤差。由于接觸式測量方法用于螺旋槳模態試驗的弊端，本文研究了利用聲信號并基于希爾伯特-黃變換(Hilbert-Huang transform，HHT)方法識別螺旋槳的模態參數。HHT包含經驗模式分解和希爾伯特變換兩步，是一種處理非線性、非平穩信號的方法[8-9]。目前，該方法在水利和建筑等行業有應用的先例。范興超等[10]使用該方法針對飛機模型的加速度響應信號進行了模態參數的識別；針對位移響應信號，魏博文等[11]基于HHT識別拱壩的固有頻率和阻尼比。在利用聲信號并基于HHT方法識別結構模態參數方面，夏茂龍等[12]研究了利用近場聲信號并基于HHT方法識別結構模態參數，并通過數值算例說明了該方法的準確性。本文基于聲學外場邊界元理論，推導了結構物在沖擊作用下的聲輻射信號與模態參數的數學模型，突破了進行模態參數識別必須采用近場聲信號的限制，并通過數值模擬進行驗證，證明了所建立數學模型的準確性，同時還驗證了該方法在較低的信噪比(signal-noise ratio，SNR)下仍具有足夠高的識別精度。在此基礎上，本文將該方法應用到螺旋槳模型的模態參數識別試驗中，成功通過該方法識別得到螺旋槳模型的前兩階模態頻率和阻尼比，明確該方法在實際工程中實用、有效。

1 理 論

1.1 結構振動聲信號與模態參數的關系

自由空間中，結構物振動聲輻射的過程可通過基爾霍夫-赫姆霍茲方程進行描述[13]，如式(1)是無空間分布聲源的情況，空間中一點的聲壓PF(r,t)可通過振動結構表面的聲壓和速度分布積分得到

(1)

式中：r′為位移矢量；nr′為結構表面法向；S為結構物表面；P(r′)為結構表面的聲壓分布；G(r,r′;ω)=exp(ik·|r-r′|)/4π|r-r′|為三維空間格林函數。對式(1)，通過邊界元方法，可表示成矩陣方程的形式，首先將場點布置于結構物表面，并作邊界元離散，得到矩陣方程式如下

(2)

(3)

(4)

結構物的振型方程可以采用有限元形式表示如下

{F(t)}

(5)

{x(t)}=[Φ]{q(t)}

(6)

(7)

(8)

[Mr]-1[Φ]T{F(t)}

(9)

式中：[Cr]=[Mr]-1[Φ]T[C][Φ]，其中[C]為阻尼矩陣；[Kr]=[Mr]-1[Φ]T[K][Φ]，[K]為剛度矩陣；[Mr]=[Φ]T[M][Φ]，[M]為質量矩陣，這些都是對角矩陣。矩陣方程式(9)中，第j階模態節點的振動方程可為

(10)

式中，{F(t)}為作用于結構物表面作用力向量，假設在k節點作用一個沖擊載荷fk(t)=δ(t)并且其他節點上沒有作用力，則式(10)可表示成

(11)

式中:ωj為第j階固有頻率；ζj為第j階模態的阻尼比；φjk為j階振型向量的第k個元素。

對于式(11)，可求解得到

(12)

由式(7)有

(13)

其中，第p個節點的速度響應可以表示為

(14)

其中

(15)

表示第j階模態振型對于節點p速度響應的貢獻。φpj,k為j階振型中節點p和k的相位差，Bpj,k=φpj/mωdj。將式(15)代入式(4)，則由于結構物第j階模態響應引起的空間中場點F的聲壓可以表示成

(16)

式(16)建立了場點F的聲信號與結構物模態參數之間的關系。

1.2 HHT方法識別模態參數

1.2.1 經驗模態分解方法

經驗模態分解方法(empirical mode decomposition, EMD)認為任何復雜的信號都是由一些不同尺度的本征模態函數(intrinsic mode function, IMF)構成。每個IMF具有以下特征：極點和零點的數量相等或者相差一個，極點之間形成的上下包絡線的平均值始終為零。由于具有這樣的性質，IMF具有良好的希爾伯特變換特性。其分解過程表達如下

(17)

式中：P(t)為原始信號；Pj(t)為從原始信號分解出來的m階IMF；r(t)為原始信號中的常量。對于式(3)所表示的聲信號，經EMD分解之后的每一個IMF并不和每一階模態響應一一對應。文獻[14]證明，先以某一階模態頻率為中心頻率對原始信號進行帶通濾波，再進行EMD，得到的第一階IMF就是結構的對應階模態響應數據。在帶通濾波處理上，文獻[15]和夏茂龍等的研究采用Chebyshev I 類帶通濾波器，本文工作中經過實踐，從減小濾波信號相位延遲考慮，采用有限脈沖響應(finite impulse response，FIR)線性濾波器。

1.2.2 模態參數識別

對采樣得到的聲壓進行希爾伯特變換及化簡過程可參考Huang等的研究和文獻[16]。最后得到聲壓解析信號的瞬時幅值和瞬時相位如下

(18)

(19)

對于式(18)等號兩邊取對數，得到

(20)

即瞬時幅值的對數和時間之間存在線性關系，斜率就是-ζjωj。并且從式(19)可見瞬時相位和時間成線性關系，斜率是ωdj。假設

-ζjωj=k1

(21)

ωdj=k2

(22)

斜率k1和k2可通過線性擬合得到，聯立

(23)

就可以計算得到第j階模態的固有頻率ωj和阻尼比ζj。

經以上分析，利用結構物的沖擊聲輻射信號并基于HHT方法進行模態參數識別的流程總結如圖1所示。

圖1 模態參數識別流程圖

2 數值模擬

在本章，通過數值模擬首先對目標槳模進行模態分析，得到槳模的前兩階固有頻率，在此基礎上，分別基于瞬態振動分析方法和瞬態聲學邊界元方法進行槳模在沖擊載荷作用下的振動響應和聲輻射仿真，最后提取一段時間的聲信號進行HHT分析，從而從數值仿真的角度評估HHT方法識別螺旋槳模態參數的精度。在槳模的瞬態振動分析中，為驗證基于聲信號的HHT算法識別模態阻尼比的準確性，人為設定前兩階模態阻尼比都為0.01。需說明的是，0.01并非是真實的阻尼比，人為設定該參數只是為了驗證算法的準確性。螺旋槳模型直徑為240 mm，盤面比約為0.74，帶有30°側斜角，材質為鋁合金，根據質材檢測報告，彈性模量為59 GPa，密度為2.7 t/m3。

2.1 槳葉的模態分析

通過有限元方法對槳葉進行正則模態計算，槳模葉根處采取剛固約束。計算結果如表1及圖2所示。可見槳葉的第一階模態振型是軸向彎曲振型，第二階模態振型是弦向扭轉振型。

表1 槳葉的前兩階固有頻率

(a) 第一階振型

2.2 槳葉在沖擊振動下的近場噪聲

在槳葉正則模態分析的基礎上，基于有限元瞬時振動響應計算方法，數值模擬槳葉中部受到單位沖擊作用時的振動響應。在完成槳葉的瞬態振動計算之后，以槳葉表面的速度分布作為聲學邊界元的邊界條件，進行瞬態聲學邊界元計算。計算總時長設置為0.05 s，計算時間步長為1/20 480 s。之后，在距離槳葉表面任意位置垂直距離2 mm處取場點1，提取聲壓信號，并通過HHT方法處理并識別模態參數。場點1采集到的聲壓信號的時域和頻域曲線如圖3所示。由于進行數值計算的總時長為0.05 s，經FFT之后的頻率間隔最小為20 Hz，從圖3中可見前兩個峰值分別位于1 260 Hz和1 800 Hz，考慮頻率分辨率，可大致確定前兩階固有頻率分別處于1 240～1 280 Hz和1 780～1 820 Hz內。

(a) 時域曲線

2.3 基于槳葉聲信號的模態參數識別

在大致確定槳葉的前兩階固有頻率的范圍之后，通過Kaiser窗函數設計法設計有限脈沖響應(FIR)濾波器。FIR濾波器具有穩定及線性相位的優點[17]。針對槳葉第一階固有頻率，設計帶通范圍為1 240 Hz～1 280 Hz，對于槳葉第二階固有頻率，設計帶通范圍為1 780～1 820 Hz，設計階數都為200。兩個濾波器的幅值曲線如圖4所示。

(a)

通過所設計的濾波器，分別對場點1的聲壓信號進行濾波處理。由于濾波后的信號存在相位延遲，還需要根據階數進行相位的修正[18-19]。場點1的聲壓信號分別經過以上兩個數字濾波器濾波及相位延遲修正之后結果如圖5所示。

分別對濾波后的信號進行經驗模態分解(EMD)，EMD分解后的第一階固有模態函數(IMF)就是對應該階模態的響應數據。分別對圖5中的兩組濾波后的信號進行EMD之后得到的第一階IMF，如圖6所示。

(a) 第一階固有頻率濾波結果

(a)

分別對圖6中的前兩階模態響應信號進行希爾伯特變換(Hilbert transform，HT)，獲得槳葉前兩階模態響應信號的瞬時幅值譜和瞬時相位譜，并對瞬時幅值譜取對數，繪制成曲線，并以最小二乘法線性擬合，結果，如圖7所示。

(a) 第一階

由圖7可見，瞬時幅值和瞬時相位曲線的端點存在嚴重的端點效應，這是由于前述數據的端點效應累計造成，包括使用FIR濾波器造成的端點數據失真和使用EMD分解造成的端點數據失真。對此，在進行最小二乘擬合時，去掉尾端部的一段數據，最終得到槳葉的前兩階模態參數如表2所示。

表2 槳葉前兩階模態頻率與阻尼比

從計算結果可見，HHT方法可以很準確地識別槳葉前兩階固有頻率和模態阻尼比，相比FFT方法，固有頻率識別準確度更高，說明該方法具有很好的頻率分辨率自適應性。此外，阻尼比識別結果非常準確，和數值解一致。

場點1位于表面附近，其聲壓信號主要受到鄰近單元的影響，在文獻[20-21]中，使用點聲源模型表示其聲壓信號，其中源強為單元體積速度。本文中，式(16)建立了空間中任意場點和結構模態參數的關系，采樣點不必要布置在近場。下面在距離槳葉較遠的位置選取編號2～6共5個場點，分別是場點1沿著軸線方向平移100 mm，200 mm，300 mm，400 mm，500 mm得到(見圖8)。以場點2～6的聲壓信號進行結構物的模態參數識別，結果如表3所示。

圖8 距離槳葉較遠的場點

表3 場點2～6聲壓信號的結構模態參數識別結果

由表3可知，基于遠場點的聲壓信號識別的槳葉模態參數依然具有非常高的準確度，其中固有頻率的誤差小于1‰，阻尼比最大誤差為1%。當然，以上數據只是基于數值模擬，所采樣的聲壓信號都是理想狀態下的，不受實際環境噪聲影響。為將該方法推廣到實際工程應用中，還需要對環境噪聲的影響進行評估。

2.4 環境噪聲對于模態參數識別結果的影響

為評估環境噪聲對于模態參數識別結果的影響，以槳葉的第一階固有頻率f1=1 253.5 Hz以及ζ1=0.01的阻尼比構建一標準衰減信號，表達式如下

(24)

以20 480 Hz的采樣頻率對該信號進行采樣，采樣時長0.05 s，得到一組標準衰減信號，另外，通過添加高斯白噪聲方式來模擬環境噪聲，并定義信噪比為1～10 240 Hz頻帶內的目標信號和噪聲信號的總能量之比，如下

(25)

式中：ST為目標信號；N為噪聲信號。通過設置白噪聲信號的幅值，可將SNR大概限制在某一數值，可定量地評估信噪比對于模態參數識別精度的影響。以下研究信噪比分別為9 dB，6 dB，3 dB，0，-3 dB，-6 dB，-9 dB，-12 dB，-15 dB情況下模態參數的識別結果。需要說明的是：由于白噪聲通過隨機數算法產生，因此，對于某個信噪比，不一定能夠精確達到，例如，目標信噪比為3 dB，實際信噪比可能是3.1 dB。由于白噪聲的偶然性，因此，對于每一個信噪比下的模態參數識別仿真，都重復100次，取結果的平均值。以上不同信噪比下模態參數識別結果如表4所示；部分信噪比的時域曲線如圖9所示。

(a) 9 dB

表4 不同信噪比的模態參數識別結果

由表4可知，隨著信噪比降低，固有頻率和模態阻尼比的識別結果精度都隨之降低。其中，模態阻尼比的誤差相對較大，相比之下，即使在較低的信噪比(-15 dB)下，固有頻率的識別結果依然可以保證比較高的精度(0.11%)。相對而言，在低信噪比下，模態阻尼比存在較大的誤差，這是因為模態阻尼比的數值比固有頻率小5個數量級，在運算過程中，同樣絕對值的誤差會對模態阻尼比的計算精度造成更大的影響。

實際工程環境中，由于無法保證測量環境的全自由場屬性以及無法避免背景噪聲，所測量得到的信號會存在一定的噪聲污染，本數值算例說明在相當低的信噪比之下，基于HHT方法依然可以較為準確地識別模態參數，證明該方法在工程應用中的實用性。

3 模型試驗

下面從試驗角度對HHT方法識別模態參數的實用性進行驗證。試驗過程中，為實現較高的信噪比，傳聲器布置于槳葉壓力面0.7倍半徑弦長中部，與表面的垂直距離約為5 mm，通過力錘敲擊槳葉梢部并同時記錄聲壓信號。試驗示意圖和實景圖如圖10所示。傳聲器采集的聲壓信號如圖11所示，采樣頻率為25 600 Hz，采樣時長為1 s。

(a)

圖11 傳聲器采集的聲壓信號

從原始信號中間位置截取一段進行分析，截取段樣本長度為2 048。對于所截取的信號，通過1.2節中羅列的步驟進行處理，最終得到槳葉第一、第二階固有頻率和模態阻尼比，如表5所示。表5中將試驗測量的固有頻率和模態數值計算結果(見表1)進行了對比，結果顯示試驗和數值結果比較接近，第一、第二階的差異分別為0.9%和2.7%。槳葉第一階模態識別數據處理過程的部分曲線圖如圖12所示，由圖12可知，對數幅值曲線存在小振幅的波動，而相位曲線則具有非常高的線性度，第二階模態識別過程的相關曲線與第一階類似。

(a) 截取的信號段

表5 螺旋槳模態參數識別結果

4 結 論

本文通過自由場環境下外部直接聲學邊界元方法的推導，建立了結構振動時外場輻射聲壓和結構物模態參數之間的關系，相比于夏茂龍等和Xu等的研究在結構物近場提取聲壓信號進行模態參數的識別分析，從理論上闡明了在遠場采集的聲信號同樣可以進行結構物模態參數的識別。在此基礎上，通過數值仿真手段驗證了所建立理論模型的遠、近場適用性以及在較低信噪比下的準確性。最后，將該方法應用在一槳模的模態參數識別試驗中，成功識別槳模的前兩階固有頻率及模態阻尼比。本文的工作表明基于HHT算法的聲信號模態參數識別方法可為槳模等精細結構的模態參數識別提供一種新的途徑。