斜下降旋翼槳渦干擾聲源定位及表面壓力

劉正江，汪文濤，林永峰，曹亞雄

中國直升機設計研究所 重點實驗室，景德鎮 333000

旋翼是直升機特有的主要升力面和飛行姿態控制舵面，同時也是直升機最大的噪聲源[1]。尤其在諸如小角度下降、拉平著陸等飛行狀態，因離地面近，其噪聲對周邊環境必然產生更強烈的擾民影響。在這些近地飛行狀態，旋翼槳-渦干擾噪聲占據直升機噪聲能量主要部分[2-3]。

針對槳-渦干擾(Blade-Vortex Interaction,BVI)噪聲形成機理、聲源位置以及傳播特性，國內外一直在持續不斷地開展著深入的研究。美國、德國和法國多家機構先后在1990年和2001年聯合發起HART I和II計劃，以B0105為試驗對象，在DNW(German-Dutch Wind)消聲風洞系統地開展了高階諧波控制(Higher Harmonic Control, HHC)、單片槳葉控制(Individual Blade Control, IBC)的降噪試驗研究[4-5]。航班直升機公司先進技術研究所為評估新設計的AT-2模型旋翼的基本參數，于2000年2月～3月在DNW開展了第2階段模型旋翼試驗，研究對比了基本型、1型和2型槳葉多種飛行狀態下HHC對槳渦干擾噪聲控制和振動控制的影響[6]。2002年4月，歐盟11個合作者共同啟動了HeliNOVI項目，相關試驗在德國航空航天研究院的DNW聲學風洞進行，主要研究目標是減小振動和降低噪聲。試驗測量得到了矩形槳尖標準槳葉的瞬時聲壓、槳葉表面壓力、動力學以及其他性能數據，研究了槳葉表面壓力和噪聲之間的關系，同時獲得了主/尾槳-渦干擾過程中的槳尖渦軌跡和槳葉的相互位置關系的定量信息，建立了一個載荷和噪聲綜合數據庫[7]。國外開展的上述項目不僅研究了槳渦干擾噪聲，還同時研究了相同試驗狀態下的槳葉載荷/表面壓力、流場以及槳葉形變。

中國如中國空氣動力研究與發展中心[8]、南京航空航天大學[9-11]以及中國直升機設計研究所[12-13]開展過旋翼噪聲試驗研究，但大多進行的是單獨噪聲試驗，同步開展相同狀態下的斜下降旋翼槳-渦干擾聲源定位試驗和槳葉表面壓力相關性試驗較少。

本文首先對槳-渦干擾狀態的槳葉表面載荷變化進行數值計算，對旋翼噪聲數據去噪方法、槳-渦干擾識別、分離以及聲源定位方法進行了研究，在直升機風洞同步開展模型旋翼斜下降槳-渦干擾聲源定位和槳葉表面壓力測量試驗，研究槳-渦干擾狀態、槳葉表面壓力載荷變化特性以及槳-渦干擾噪聲源特性，從時間和空間維度上分析槳-渦干擾噪聲和槳葉表面壓力載荷之間的相關性。

1 槳渦干擾狀態槳葉表面載荷數值計算

1.1 自由尾跡計算槳葉表面載荷方法

采用自由尾跡方法進行斜下降槳渦干擾狀態的槳葉表面氣動載荷計算[14]。自由尾跡方法將旋翼尾跡離散為渦元，而渦元可以隨當地氣流運動，能很好地描述旋翼運動尾跡。計算中將尾跡分為近尾跡和遠尾跡分別處理，由近尾跡過渡至遠尾跡的渦齡角為60°。近尾跡區域內槳尖渦處于卷起過程中，尾跡保持渦面的形式；遠尾跡槳尖渦完全形成，以線渦單元的形式對槳尖渦進行處理，并采用渦核模型消除槳尖渦附近計算奇異性。計算時先將槳葉離散為有限的氣動分析單元，再采用Biot-Savart公式計算槳葉在不同展向位置和不同方位角的誘導速度，進而實現不同剖面槳葉的非定常載荷計算。Biot-Savart公式為

(1)

式中：dv為微段誘導速度；Γ為微段渦量；dσ為槳尖渦渦線微段；r為微段到周圍流場任意點的距離的矢量；r為微段到周圍流場任意點的距離。

1.2 自由尾跡計算結果

風洞試驗采用4 m直徑量級測壓模型槳葉，槳轂為四支臂鉸接式，槳葉半徑R為2 100 mm，槳葉弦長D為140 mm，槳尖形狀為矩形，采用OA213翼型。計算時選取槳葉表面載荷一個典型觀測剖面位置：槳葉徑向0.85R處。同時，旋翼轉速為749 r/min，前進比為0.191(對應風速為30 m/s)。數值計算的槳葉和狀態與風洞試驗保持一致，槳葉表面氣動載荷計算狀態如表1所示。槳葉0.85R剖面的氣動載荷計算結果如圖1所示。

表1 計算參數Table 1 Computation parameters

圖1 計算槳渦干擾狀態槳葉表面載荷Fig.1 Computed rotor surface load in blade-vortex interaction

從圖1可以看出，槳葉表面載荷在320°方位角附近(圖1中虛線圓)振蕩劇烈，表明槳葉在此方位角(后行側)附近發生了較強烈的槳-渦干擾現象。槳葉表面載荷在后行側320°方位角附近會產生槳-渦干擾，這與國外采用相近似參數試驗件及狀態得出的結果基本一致[5]。

2 槳-渦干擾噪聲數據處理方法

2.1 改進型整周期時域同步平均去噪方法

旋翼旋轉產生的氣動噪聲具有典型的以轉速為基頻的周期性特征，因此數據處理時首先對采集到的旋翼噪聲信號進行整周期平均以剔除非旋轉噪聲干擾成分。整周期平均的主要思路為：利用同步采集的0方位脈沖鍵相信號作為表面壓力和槳-渦干擾噪聲數據塊的起始點和結束點，將連續采集到的數據分割成整周期信號，然后將截取的整周期信號進行同序點平均。該平均方法充分利用了信號的相關性，不僅可以有效抑制隨機性干擾噪聲并提高旋翼噪聲的信噪比，同時能較好地保持旋翼噪聲信號幅值和相位特征信息[15-17]。

設測量信號由信號成分S(t)(確定性分量)和噪聲干擾n(t)組成，則進行N次測量后，在某特定時刻ti經同步平均后的信噪比為

(2)

式中：Sj和nj分別為第j次采集時候的信號成分和干擾噪聲信號成分。

試驗時測量采用定頻采樣，采樣率為每秒32K個采樣點。高采樣時轉速對應到單個周期后會有輕微波動，體現在每圈數據長度會稍有差別，如果不進行相應處理而是直接進行整周期同序點平均，會造成相位錯亂并最終導致后續聲源定位出現較大誤判。針對這種情況，對算法進行了如下改進：先按照整周期時域同步平均的方法將采集的若干秒數據截取成多個整圈數據，再統計出每圈數據的長度，得到數據長度多數值，最后從中挑選出長度為多數值的若干圈數據(一般不少于50 r，否則需增加采集數據時間)進行整周期同序點的同步平均。

2.2 基于小波變換槳渦干擾信號識別和分離方法

槳-渦干擾狀態下的槳葉表面壓力和槳-渦干擾噪聲信號具有典型的脈沖特性，采用常規的濾波方法無法保證時域波形和相位不失真，而失真必將導致聲源定位計算結果嚴重失實。基于小波變換的小波包處理不僅具有良好的時頻局部細化能力，而且不會造成原有信號特征丟失。小波包分解是將頻帶進行多層次劃分，再通過多分辨率分析方法對沒有細分的頻段進行進一步分解，從而提高時頻分辨率。而根據熵理論及意義，不同類型的旋翼噪聲在各個頻段的能量是不一樣的，熵值隨之也不同，熵值對于信號的時變比較敏感，由其定義的小波包特征熵可表征信號復雜度在時域的變化情況以及豐富的頻域特征[18-20]，因此選用其對槳-渦干擾噪聲進行識別和分離處理。

利用小波包進行槳-渦干擾信號識別的具體步驟如下：

第1步信號分解

在LabVIEW環境下利用小波包工具箱提供的函數對整周期平均去噪后的槳-渦干擾信號進行n層小波包分解，分別提取第m(m≤n)層從低頻到高頻的k個頻帶信號；

第2步信號重構

對第1步之后得到的若干頻帶的小波包系數xj,k進行重構，得到k個小波包重構信號Sm,k(t)，其中各個重構信號分別包含了信號從低頻到高頻的各個頻帶的信息，信息量既無冗余，也不疏漏。對于m層的小波包分解而言，它是按照從最低頻到最高頻等頻率間隔將信號分解到k個頻帶上。這里假設k=8，m=3，n=10。

第3步計算各頻帶信號的能量

設S3,k(t)(k=0,1,…,7)對應能量E3,k(k=0,1,…,7)，則有：

(3)

式中：xj,k(j=1,2,…,7;k=0,1,…,n)表示重構信號S3,k(t)離散點的幅值。對各個節點能量做處理組成向量，設節點向量為

E=[E3,0,E3,1,E3,2,E3,3,E3,4,E3,5,E3,6,E3,7]

(4)

再將各個節點能量做處理組成新的向量：

C=[C3,0,C3,1,C3,2,C3,3,C3,4,C3,5,C3,6,C3,7]

(5)

第4步構造信號小波包特征熵向量

(6)

式中：H3,k表示經過3層小波包分解后在8個頻帶的小波包特征熵。分別求得8個小波包特征熵，進而以8個小波包特征熵為元素可以構成一個特征向量表示為

H=[H3,0,H3,1,H3,2,H3,3,H3,4,

H3,5,H3,6,H3,7]

(7)

2.3 波束成形聲源定位方法及軟件實現

常用的聲源定位方法主要有聲強法、近場聲全息(Nearfield Acoustical Holography，NAH)和遠場波束形成(BeamForming)等，其中波束形成聲源定位方法因測點需求少、計算速度快、中高頻分辨率好、測量距離遠、覆蓋面積大以及易于布置等特點獲得廣泛工程應用[21-22]。

波束形成算法通常分兩步來計算得到聲源位置：首先計算不同麥克風對之間的時延；然后根據計算得到的時延值以及麥克風陣列的幾何形狀，采用幾何法或插值法等方法來確定聲源具體的空間位置。

2.3.1 bartlett時延算法

選用bartlett法計算時延，bartlett法是一種基于平均周期圖法的功率譜估計方法，基本原理是先將N點有限長度序列分段后求周期圖并進行平均，計算平均數據段的功率譜，再將功率譜進行多項式擬合得到擬合后的多項式系數矩陣，最后計算矩陣特征值即獲得最后所需的時延值[23]。具體步驟如下：

首先，記輸入信號數據為Y(n)，取q段平均，每段數據長度為L，第i段數據記為

Y(i)(n)=Y(n+(i-1)L)

(8)

式中：i=1,2,…,q;n=1,2,…,L。

第i段的周期圖為

(9)

則平均后的功率譜為

P(k)=10lgPxx(k)

(10)

式中：Pxx為自功率譜函數。

然后對其進行多項式擬合，即：

(11)

式中：q0為多項式的項數；αq為常數項。則功率譜特征可表示為A=[α0,α1,…,αq0]。

再利用式(12)計算得到矩陣T：

T=ATA

(12)

最后計算矩陣T的特征值即可得到時延di0。

2.3.2 基于球面插值法的聲源定位方法

在計算得到陣列傳聲器時延di0的基礎上，后續確定聲源位置的方法大體上分兩種，一種是二維空間搜索方式，另一種是幾何計算方法。二維空間搜索方法雖然可以精確定位空間聲源位置，但當所用陣列傳聲器較多時，其搜索非常消耗計算時間和資源，因此在工程應用上受計算條件限制一般只能用于事后處理。而為提高計算效率以保證試驗時聲源定位結果展示具有一定的實時性，對比分析后選用了幾何方法中綜合效能最優的球面插值法進行聲源定位計算[24-26]。

圖2為球面插值法的空間定義。根據矢量幾何和三角形三邊關系，有：

新泰地處魯中腹地，全市有林地面積達到117.1萬畝，森林覆蓋率達到39.8%，是山東省森林資源比較豐富的縣之一。近年來，新泰市以明晰所有權、穩定承包權、放活經營權為核心，吸引社會資本進山入林，培育林業新型經營主體，打造林地集體所有、家庭承包、多元經營新格局。新泰被評為國家園林城市、全國綠化模范市、全國經濟林產業區域品牌建設試點單位，“5+5”造林和經營模式在全省推廣。

(13)

(14)

式中：RS為聲源到傳感器0的距離；Ri為傳感器i到傳感器0的距離；ri為聲傳感器i的空間向量ri=[xi,yi,zi]；rS為聲源的空間向量。

由于di0是由bartlett時延計算得到的，與真實值有一定的誤差，為此引入誤差函數：

(15)

將式(15)得到的N個方程寫成誤差矩陣為

ε=δ-2RSd-2SrS

(16)

式中：

R*RT-d*dT

(17)

d=[d10,d20,…,dN0]T

(18)

(19)

圖2 球面插值法空間定義Fig.2 Space definition of spherical interpolation

其中：R為由Ri組成的N× 1矩陣；d成為由di0為組成的N× 1矩陣。

式(16)可表示成：

(20)

式中：

(21)

采用最小平方誤差估計，即使誤差的平方和最小，表示為

(22)

為了使eLS達到最小，需要

(23)

寫成矩陣分塊形式為

(24)

簡化為

(25)

式中：

(26)

式中：I為單位矩陣。

定義投影矩陣：

(27)

則聲源位置的估計值可以化簡為

rS,SI=(STPdS)-1STPdb

(28)

2.3.3 槳渦干擾噪聲分析軟件實現

軟件采用LabVIEW、高級信號處理工具包和聲音振動工具包開發，支持的前端采集硬件包括德國Bustec公司的3424系列VXI總線動態采集卡和美國NI公司的4472系列PXI/PCI總線動態采集卡。系統運行后可以自動識別在線硬件并生成通道參數配置表，用戶可通過該表對傳感器系數、三維空間坐標等參數進行配置(見圖3)。主界面(見圖4)中可以選擇時間歷程、幅值譜、倍頻程以及聲源定位云圖等多種顯示方式，用戶可修改角度范圍、比例系數等聲源定位計算參數，聲源定位云圖中給出最大聲源坐標位置、最大聲源頻率和總聲壓級等信息。

圖3 通道配置界面Fig.3 Channel configuration interface

圖4 聲源定位軟件主界面Fig.4 Noise source location software main interface

3 槳葉表面壓力傳感器和聲陣列布置

3.1 槳葉表面壓力傳感器布置

試驗前，在4號槳葉上表面弦向0.1D，徑向0.90R、0.85R和0.80R3個位置分別布置KULITE的LQ-062微型表面壓力傳感器(見圖5和圖6)。槳葉加工時按照壓力傳感器大小預留圓形凹槽，壓力傳感器采用CN(Cyanoacrylate)膠固化在圓形凹槽內，壓力傳感器線纜走線采用和槳葉表面粘貼應變花相同的走線工藝，壓力傳感器的恒溫電阻以上下溝槽壓緊方式固定在槳轂中心處的輕質鋁盒中。

圖5 表面壓力測點布置Fig.5 Arrangement of surface pressure measuring points

圖6 表面壓力測點實物圖Fig.6 Photograph of surface pressure measuring points

3.2 麥克風聲陣列選擇及布置

常用的麥克風陣列包括十字陣列、矩形陣列、同心圓陣列以及螺旋陣列等，相對來說螺旋陣列具有最優的空間分辨率(即定位精度)，因此選擇螺旋陣列，陣列有效直徑為3.0 m。

受到試驗場地空間位置以及陣列直徑限制，同時采用的是遠場聲源定位算法，綜合考慮上述因素后決定將聲陣列布置在風洞流場外槳盤平面正下方偏前位置處，參考相關資料[9]陣列放置在該位置只能進行斜下降狀態后行側的槳-渦干擾噪聲聲源定位。為盡量降低地面聲反射影響，同時又便于試驗時的陣列和傳聲的安裝以及現場標定位置調整，專門設計了一個高度為1.8 m的可移動陣列安裝支架(見圖7)。該安裝支架可以通過滑輪推動，支架上部聲陣列安裝平面可以繞軸旋轉已調節聲陣列的俯仰角。試驗前根據設定的試驗狀態，取試驗時斜下降狀態主軸傾角(此時，相當于槳盤平面后倒，見圖8)中間值5.5°，將聲陣列前移，陣列平面前倒到5.5°，并使旋翼主軸和陣列平面中心垂直軸在一條直線上，從而保證陣列平面和試驗時的槳盤平面盡量(不含周期變距)平行。試驗前安裝聲陣列支架以及地面包裹、鋪設吸聲海綿。

圖7 聲陣列實物圖Fig.7 Photograph of microphone array

圖8 聲陣列風洞布置示意圖(左視)Fig.8 Schematic diagram of microphone array in wind tunnel (left view)

聲陣列布置空間位置示意圖見圖8。5 500 mm為槳轂中心到駐室進氣段最近點水平距離，實線聲陣列擺放的位置為槳盤平面正下方，虛點表示現場標定時其位置朝進氣段前移998 mm，即移至槳盤平面正下方偏前位置。

4 聲陣列圖像畸變修正及現場校準

4.1 圖像畸變修正算法

由于3 m直徑螺旋聲陣列不是完全剛性平面，無法精確保證所有傳聲器都在一個平面內，即螺旋聲陣列平面和槳盤平面實際上存在一個很小的夾角，致使聲源定位圖像有一定的畸變失真，因此需要對其進行聲源定位圖像畸變矯正。

假設P(u,v)和Pd(ud,vd)分別是失真和無失真的圖像坐標，兩者之間的坐標位置關系如下：

(29)

式中：(u0,v0)為坐標原點(聲源定位云圖中心點)；Ki為徑向畸變系數；rp為圖像點到圓心的距離。實際使用只需考慮前兩個畸變項K1和K2。

4.2 聲陣列現場校準及精度分析

試驗前，將旋翼槳盤平面后倒5.5°，通過在0°、90°和270°槳尖位置附近的定頻全向喇叭(1 000 Hz 單頻聲)對聲陣列參數進行了校準和圖像畸變修正，最終結果見圖9。

圖9 聲陣列參數校正后聲源定位結果云圖Fig.9 Cloud image of sound source localization results after acoustic array parameter correction

從表2可以看出，經過校準和圖像畸變修正后，聲源定位最大方位差為2.883°、最大徑向位置差為0.010R，滿足試驗對上述值誤差要求。

表2 現場校準聲源定位誤差分析Table 2 Error analysis of noise source location calibration

5 典型試驗數據分析

5.1 槳渦干擾狀態BVI噪聲數據分析

選取一組典型試驗狀態數據進行分析，試驗狀態為旋翼轉速749 r/min、風速30 m/s、主軸傾角5.5°，取狀態穩定后的10 s數據，按2.1和2.2節介紹的數據處理方法，得到平均后的一圈數據，再進行BVI噪聲的小波識別和分離，處理前后的結果見圖10。

對比圖10(a)和圖10(b)可以得到，采用2.1和2.2節介紹的方法能夠清晰識別并提取出槳渦干擾噪聲，同時能夠較好地保留槳渦干擾噪聲的幅值和相位波形特征。

圖10 聲陣列噪聲整圈聲壓歷程Fig.10 One rotating period noise data of microphone array

5.2 槳渦干擾狀態槳葉表面壓力數據分析

選用和5.1節相同試驗狀態數據以及數據處理方法，對旋翼轉速749 r/min、風速30 m/s的不同主軸傾角試驗狀態(見表3)的槳葉表面壓力數據進行平均及識別和提取處理。同時，已知粘貼微型表面壓力傳感器的4號槳葉方位角為266°，轉速0脈沖傳感器的方位角為356°，據此對提取出的表面壓力數據進行方位角修正。圖11給出典型試驗狀態下0.90R、0.85R和0.80R剖面的一整圈(0°～360°方位)的時間歷程曲線，軌跡角為6.5°。

從圖11以及表3可以看出：

1) 槳渦干擾引起的槳葉表面壓力變化主要集中在320°方位和60°～70°方位附近，前者為后行側槳渦干擾引起的槳葉表面壓力變化，后者為前行側槳渦干擾引起的變化。

表3 槳葉表面壓力數據Table 3 Rotor surface pressure data

2) 從波形特性對比來說，后行側的峰值普遍明顯比前行側大，后行側表面壓力時域波形具有很明顯的槳渦干擾脈沖特性，而前行側槳渦干擾脈沖特性相對來說不是非常突出。

3) 大致上最大表面壓力峰值出現在主軸傾角5.5°、后行側320°方位的0.80R剖面位置處；同時隨主軸傾角繼續增加，最大表面壓力峰值位置逐漸偏向0.90R剖面。

5.3 槳渦干擾狀態BVI聲源定位數據分析

采用2.3節介紹的BeamForming波束成形法，對5.1節提取出的BVI噪聲數據進行聲源定位處理，得到最終聲源定位云圖。圖12分別給出了主軸傾角4.5°、5.5°和6.5°聲源定位云圖，表4給出各主軸傾角下的最大聲源位置、總聲壓級及峰值頻率等信息。

圖11 表面壓力時間歷程曲線Fig.11 Time history curves of surface pressure

圖12 聲源定位結果云圖Fig.12 Noise source location maps

表4 BVI聲源定位數據Table 4 BVI noise location data

從圖12及表4中可以看出：

1) 試驗時聲陣列放置在旋翼正下方偏前側，定位出的聲源位置主要在槳葉后行側，理論上放置在該位置的陣列傳聲器也只能拾取到后行側槳渦干擾噪聲(后行側槳渦噪聲主要向正下方傳播，前行側槳渦干擾噪聲主要向前下方傳播)，相關參考資料[7]已驗證了該結論。

2) 各試驗狀態下定位出的最大聲源位置主要分布在310°～330°方位角區間，而且多數集中在320°方位角附近，這和5.2節槳葉表面壓力測量分析得到的結果以及1.1節數值仿真計算結果相符。

3) 各試驗狀態的最大聲源位置總聲壓級中，以主軸傾角5.5°時為最大，此時的槳葉徑向位置大致在0.81R左右，這和5.2節槳葉表面壓力測量分析結果0.80R基本一致。

4) 隨主軸傾角進一步增加(5.5°以后)，最大聲源位置由0.80R逐漸偏向0.90R，這和5.2節槳葉表面壓力測量分析結果也一致。

6 結 論

1) 針對斜下降狀態槳渦干擾噪聲具有的整周期特性，首先設計了一套改進整周期時域同步平均去噪方法去除和槳葉旋轉倍頻不相干的干擾噪聲，然后在此基礎上采用多層小波包分解和重構算法，從去噪后的旋翼噪聲中識別并提取出槳渦干擾噪聲，很好地保留了槳渦干擾噪聲的幅值和相位等波形特征。

2) 采用bartlett時延算法和球面插值算法對提取出的槳渦干擾噪聲進行聲源定位分析，獲取聲源定位云圖及最大聲源位置坐標、頻率和聲壓級等信息。

3) 設計并開展了模型旋翼風洞斜下降槳渦干擾聲源定位和表面壓力相關性試驗，試驗分析結果表明槳盤下方陣列只能定位出后行側的槳渦干擾聲源位置，且隨主軸傾角增大，槳渦干擾聲源位置逐漸向槳尖靠近。

4) 試驗數據表明后行側槳渦干擾表面壓力峰值明顯比前行側大，這和數值計算結果基本一致；而后行側相對較大初步分析是槳尖速度和前進比較小以及槳葉翼型較厚等因素綜合造成的。

5) 從時間和空間上對槳渦干擾表面壓力和最大聲源位置進行了相關性分析，兩者主要結論相吻合，并同時表明小速度斜下降狀態槳渦干擾發生的位置主要集中在320°方位左右，從而驗證了所實現的算法以及試驗方法的合理性和實用性。