朱旭程，馬 強，侯志強

(海軍航空工程學(xué)院 飛行器工程系，山東 煙臺 264001)

旋翼是直升機飛行的關(guān)鍵動力部件，在槳葉周期變距和氣動交變載荷作用下，槳葉、槳轂、減擺器及變距拉桿等會出現(xiàn)磨損與疲勞失效故障。旋翼故障為高等級危險故障。準(zhǔn)確監(jiān)測旋翼的損傷狀況，對提高直升機飛行安全極具現(xiàn)實意義。Azzam等［1］研究旋翼損傷監(jiān)測，提出通過槳葉相對位移與機體振動識別旋翼故障方案。Ganguli等［2-3］分析了典型損傷對旋翼氣彈響應(yīng)信號影響，運用神經(jīng)網(wǎng)絡(luò)、模糊推理等人工智能技術(shù)建立旋翼故障診斷算法。在飛行條件下對旋翼損傷進行檢測尚存技術(shù)瓶頸［4］:① 在高速旋轉(zhuǎn)的旋翼上直接加裝各種數(shù)據(jù)采集系統(tǒng)遭技術(shù)限制。直升機健康及使用監(jiān)測系統(tǒng)(Health and Usage Monitoring System，HUMS)主要通過在機身上安裝傳感器獲得旋翼狀態(tài)數(shù)據(jù)，可用信號種類較有限;②因槳距時變、氣彈耦合及非定常尾流等因素影響，用傳統(tǒng)數(shù)據(jù)分析方法難以獲得對損傷敏感而對正常狀態(tài)擾動、噪音、建模誤差不敏感的旋翼損傷觀測信號，缺乏有效的旋翼損傷跟蹤方法。為此，本文用源自混沌理論的相空間重構(gòu)技術(shù)，將旋翼狀態(tài)監(jiān)測數(shù)據(jù)變換到更高維數(shù)的狀態(tài)空間內(nèi)，以不同時間尺度研究旋翼系統(tǒng)的損傷演化過程。以旋翼實測相軌跡與基準(zhǔn)旋翼狀態(tài)短時預(yù)測模型計算相軌跡之差異作為新?lián)p傷觀測信號，提出新的旋翼損傷跟蹤方法。此法可根據(jù)單個旋翼狀態(tài)監(jiān)測信號的時序數(shù)據(jù)在相空間中重建旋翼系統(tǒng)的動力學(xué)本質(zhì)，對旋翼的狀態(tài)監(jiān)測更有效。

1 旋翼在重構(gòu)相空間的損傷演化特性

1.1 旋翼氣彈有限元模型

為對比分析健康與損傷旋翼的動態(tài)特性，研究損傷對槳葉位移、槳轂力等旋翼狀態(tài)監(jiān)測信號影響，建立旋翼氣彈有限元模型。由Hamilton能量變分原理得:

式中:Π為系統(tǒng)總能，T為動能，U為應(yīng)變能，W為外力虛功，氣動載荷采用葉素理論建立，旋翼誘導(dǎo)速度采用White-Blake入流模型［5］計算。用有限元方法離散能量積分，用分部積分法建立旋翼氣彈有限元方程為:

式中:M、C、K、F分別為慣性矩陣、阻尼矩陣、剛度矩陣及載荷向量，在定常平飛狀態(tài)下為方位角的周期函數(shù)。旋翼氣彈有限元方程可采用時域有限元方法化為代數(shù)方程求解，或?qū)⑵滢D(zhuǎn)換為非線性狀態(tài)方程用數(shù)值積分方法求解。

旋翼故障模式種類較多，不同類型的故障對系統(tǒng)參數(shù)影響不同。采用復(fù)合材料槳葉剛度疲勞損傷演化模型為:

圖1 槳葉剛度損傷變量擴展曲線Fig.1 Damage accumulation curve of blade stiffness

直升機HUMS系統(tǒng)能直接提供的旋翼狀態(tài)監(jiān)測信號主要有二類:一類為槳葉揮舞、擺動及扭轉(zhuǎn)位移，可通過在機身上安裝激光多普勒測振儀實時監(jiān)測［6］;另一類為槳轂力及力矩，可通過機身上的應(yīng)變計或振動傳感器結(jié)合飛行參數(shù)獲得。

以槳尖位移、槳轂力及力矩作為旋翼典型狀態(tài)監(jiān)測信號，將旋翼氣彈有限元模型與槳葉損傷演化模型相結(jié)合，可計算分析槳葉損傷對槳尖位移及槳轂載荷影響。本文以BO-105直升機旋翼為參考設(shè)置實例旋翼基本參數(shù)為:槳葉數(shù)4，旋翼半徑R=4.94 m，旋翼轉(zhuǎn)速 Ω=40 rad/s，槳葉密度m=6.46 kg/m，槳葉弦長0.055 R，弦向剛度 0.026 8 mΩ2R4，揮向剛度 0.010 8 mΩ2R4，扭轉(zhuǎn)剛度0.006 15 mΩ2R4，洛克數(shù)5.2，槳葉升力線斜率6.0，型阻系數(shù) 0.002，誘導(dǎo)阻力系數(shù) 0.2，俯仰力矩系數(shù)0.00。飛機在前進比0.3、旋翼拉力系數(shù)0.004 9條件下定常平飛。旋翼仿真程序用Matlab軟件建立(圖2)，主程序先調(diào)用公式推理模塊執(zhí)行符號積分及微分，算出各單元慣性矩陣、剛度矩陣及載荷向量，槳葉揮舞、擺動位移形函數(shù)用3次Hermitian多項式，扭轉(zhuǎn)位移形函數(shù)用2次Lagrangian多項式以保證彎、扭力矩變分的線性一致。為加快計算收斂過程引入比例阻尼，載荷向量含環(huán)量、非環(huán)量氣動載荷及節(jié)點位移的非線性項，誘導(dǎo)速度用內(nèi)嵌的迭代過程計算。主程序在龍格-庫塔積分步中循環(huán)調(diào)用右側(cè)函數(shù)模塊求解狀態(tài)方程積分。右側(cè)函數(shù)模塊消除單元矩陣及向量中的槳葉方位角等符號變量后，裝配生成整體矩陣并用高斯消元法求出狀態(tài)增量，最后采用求和法計算出槳轂載荷。槳葉在健康(t/tF=0)與損傷(t/tF=1/3)狀態(tài)下仿真產(chǎn)生的槳尖位移、槳根力及力矩如圖3所示。圖中，w為槳尖揮舞位移，v為槳尖擺動位移，F(xiàn)z為槳轂力，Mx為槳轂力矩，n為樣本長度。

圖2 旋翼仿真程序流程Fig.2 Flowchart of rotor simulation procedure

1.2 旋翼狀態(tài)監(jiān)測信號相空間重構(gòu)

旋翼損傷監(jiān)測的基本問題為能找到對損傷敏感而對系統(tǒng)狀態(tài)變化、干擾噪聲不敏感的觀測信號。而槳葉損傷與氣彈耦合現(xiàn)象會引起氣動載荷沿槳葉結(jié)構(gòu)分布發(fā)生變化［7］，隨著槳葉損傷程度的擴展，在揮舞、擺動、扭轉(zhuǎn)方向的固有頻率互相穿越與交叉，使槳尖位移與槳轂載荷峰值等特征信號在單調(diào)性關(guān)系上出現(xiàn)突變，難以據(jù)該特征信號辨識出損傷的變化趨勢。源自混沌理論的相空間重構(gòu)技術(shù)可在不同維數(shù)空間內(nèi)、以不同時間尺度研究演化系統(tǒng)的動力學(xué)特性。只要時間尺度與空間維數(shù)選擇合適，即可望通過相軌跡的移動變化定量跟蹤系統(tǒng)損傷的演化過程。將某監(jiān)測信號時序數(shù)據(jù){y(i):i=1，2，…，M}進行時間坐標(biāo)延遲，構(gòu)造一狀態(tài)點列 {x(n):n=1，2，…，N}。其中，x(n)= ［y1(n)，y2(n)，…，ym(n)］T，yk(n)=y［n+(k-1)τ］，m為嵌入維數(shù)，τ為延遲時間。該點列描述系統(tǒng)在m維狀態(tài)空間中的相軌跡，由其構(gòu)成的空間稱為相空間，此過程稱相空間重構(gòu)。據(jù)嵌入理論［8］，若重構(gòu)系統(tǒng)與原動力學(xué)系統(tǒng)的動態(tài)響應(yīng)可限制在相同維數(shù)的吸引子流形上，即使重構(gòu)系統(tǒng)與原系統(tǒng)維數(shù)相差較大，也可在拓撲關(guān)系或微分意義上屬同一等價(微分同胚)，具有相同的動力學(xué)性質(zhì)與幾何性質(zhì)，系統(tǒng)的許多特性如分維、K熵、李指數(shù)等均相同。

為直觀表明旋翼在相空間中的損傷演化過程，將帶槳葉損傷BO-105旋翼的槳尖位移及槳轂載荷數(shù)據(jù)進行相空間重構(gòu)。相空間重構(gòu)時，延遲時間與嵌入維數(shù)是兩主要參數(shù)，采用平均互信息法及假近鄰法確定［9］，計算流程見圖4。求得延遲時間為30，嵌入維數(shù)為4，信號相軌跡在三維空間中的投影見圖5。由圖可直觀看出旋翼系統(tǒng)的穩(wěn)態(tài)相軌跡(吸引子)在重構(gòu)空間中被完全展開，通過相軌跡隨槳葉損傷擴展而出現(xiàn)的漂移及變形等幾何變化可直觀了解與跟蹤旋翼損傷演化過程。

圖5 相軌跡在三維空間中的投影Fig.5 3D views of reconstructed phase space trajectories

2 重構(gòu)相空間跟蹤旋翼損傷方法

基于相空間狀態(tài)預(yù)測殘差的旋翼損傷跟蹤程序(圖6)主要包括損傷觀測信號生成模塊與損傷狀態(tài)優(yōu)化估計模塊。前者先據(jù)健康旋翼在相空間內(nèi)的運動軌跡采用Volterra預(yù)測器［10］建立基準(zhǔn)旋翼短時狀態(tài)預(yù)測模型，系統(tǒng)運行時用現(xiàn)場采集數(shù)據(jù)及Volterra預(yù)測器估計基準(zhǔn)狀態(tài)，與實測狀態(tài)之間的差異即短時狀態(tài)預(yù)測殘差，對一條記錄內(nèi)的殘差進行統(tǒng)計、組合可形成殘差特征向量。而后者將生成的殘差特征向量代入觀測方程得到槳葉損傷觀測值，用槳葉損傷演化方程計算Unscented變換［11］，產(chǎn)生西格瑪點及其權(quán)重與狀態(tài)轉(zhuǎn)移點，調(diào)用預(yù)測步計算損傷變量均值及方差的一步估計;調(diào)用修正步估計觀測變量均值、方差，計算增益矩陣并修正損傷變量的均值估計與方差。

圖6 旋翼損傷跟蹤程序原理圖Fig.6 Rotor damage tracking system framework

2.1 損傷觀測信號生成

在旋翼系統(tǒng)狀態(tài)空間中，可用損傷旋翼與健康旋翼軌跡之間的短時偏差作為系統(tǒng)損傷觀測信號。設(shè)旋翼在某任意時刻狀態(tài)為xn，損傷為φ，經(jīng)過短時間τ后旋翼狀態(tài)變?yōu)閤(xn，τ，φ)。若按健康系統(tǒng)(φ=φR)估計其狀態(tài)則應(yīng)為x(xn，τ，φR)，兩狀態(tài)之間的殘差為:

它在xn，τ一定條件下對損傷具有近似線性可觀性。

其中:

用式(5)計算出短時狀態(tài)殘差后，即可研究靈敏度S(xn，τ)隨系統(tǒng)初始狀態(tài)及損傷大小的變化情況，具體分析高階小量影響。

分析實例旋翼槳尖揮舞位移的短時殘差靈敏度隨初始狀態(tài)和槳葉剛度損傷的變化情況。在健康旋翼的狀態(tài)空間中隨機選擇16個點作為計算初始狀態(tài)，在槳葉壽命周期(0≤t/tF＜1)內(nèi)等間隔取60個時間點，用式(5)計算出旋翼槳葉揮舞位移的短時狀態(tài)殘差隨使用壽命和槳葉損傷大小的變化曲線見圖7。左圖中16條槳尖揮舞位移的短時殘差ew的曲線對應(yīng)于16個初始狀態(tài)，在初始狀態(tài)一定的情況下，槳尖位移的短時殘差是損傷變量的光滑連續(xù)函數(shù)，且與槳葉損傷演化曲線非常一致。右圖表明，短時殘差靈敏度只與系統(tǒng)的初始狀態(tài)有關(guān)，而不隨損傷的大小有明顯變化，如果對多個狀態(tài)點的短時殘差進行平均，則可以得到一種相對于系統(tǒng)狀態(tài)和損傷都有比較穩(wěn)定靈敏度的損傷觀測信號。

圖7 短時狀態(tài)殘差與壽命(左)和損傷(右)的關(guān)系Fig.7 Curves of short time state error vs.life time(left)and damage variable(right)

表面上，有了旋翼氣彈有限元模型即可計算x(xn，τ，φR)，但該方法在計算短時狀態(tài)殘差時要求能直接測量旋翼氣彈有限元模型中的全部狀態(tài)量，需實時采集許多節(jié)點處的槳葉位移數(shù)據(jù)，這在實際應(yīng)用中較難進行?？紤]到原物理系統(tǒng)在一定意義下與相空間中的重構(gòu)系統(tǒng)存在某種拓撲等價關(guān)系，則可在基準(zhǔn)相空間中利用基準(zhǔn)狀態(tài)預(yù)測模型生成新的狀態(tài)殘差ex，作為物理空間中狀態(tài)殘差e^x的替代，即，據(jù)健康旋翼在相空間內(nèi)的運動軌跡，用非線性Volterra級數(shù)建立基準(zhǔn)旋翼短時狀態(tài)預(yù)測模型。監(jiān)測時將旋翼狀態(tài)采集數(shù)據(jù)運用延遲嵌入方法重構(gòu)出實測相軌跡，并采用Volterra預(yù)測模型算出參考軌跡，并將其與實測相軌跡之間的差異統(tǒng)計平均生成損傷觀測信號。由于傳感器測量誤差、Volterra預(yù)測模型在不同區(qū)域內(nèi)的精度變化、靈敏度在相軌跡上的起伏誤差等隨機干擾影響，單個短時狀態(tài)殘差具有較大不確定性。為提高損傷檢測的魯棒性，用一條記錄內(nèi)所有殘差的平均值ek作為損傷觀測特征向量?；赩olterra預(yù)測模型的損傷觀測信號生成方法為:

其中:y為傳感器輸出序列;Pha為相空間重構(gòu)算子;m，τ分別為嵌入維數(shù)及時延參數(shù);x為相點;xR為預(yù)測器輸出參考點;ex為狀態(tài)預(yù)測殘差;ek為殘差均值(損傷觀測信號);Nk為記錄長度;Vol為Volterra預(yù)測器;p為預(yù)測器階數(shù)。P階截斷及形式的Volterra預(yù)測器可表示為:

式中:

其中:hp(i1，…，ip)為p階Volterra核;X(n)為預(yù)測器輸入信號矩陣;w為預(yù)測器系數(shù)向量，可據(jù)健康旋翼的相軌跡用NLMS自適應(yīng)算法［11］確定:

式中:μ為收斂控制參數(shù)(0＜μ＜2)。

圖8 損傷觀測信號計算流程Fig.8 Flow chart of damage measurement prediction

為說明采用Volterra預(yù)測模型生成損傷觀測信號方法的有效性，用損傷觀測信號生成程序(圖8)對實例旋翼在不同損傷狀態(tài)下的槳葉位移及槳轂載荷進行分析。在槳葉壽命周期內(nèi)(0≤t/tF＜1)，等時間間隔采集60個損傷程度(0≤φ＜1)下的仿真記錄，計算得到各記錄殘差均值隨使用壽命的變化如圖9所示。與圖1對比可見，各殘差均值分量與槳葉損傷變量的變化趨勢較一致，說明殘差主要來自槳葉損傷的變化而非模型計算誤差。對比分析殘差均值和槳葉損傷之間的對應(yīng)關(guān)系，容易驗證兩者有較顯著的線性關(guān)系。

圖9 殘差均值向量隨時間的變化Fig.9 Curves of residual tracking metric

2.2 槳葉損傷估計

旋翼系統(tǒng)偏離正常狀態(tài)時，殘差均值ek可近似線性反映出槳葉的損傷程度。設(shè)由殘差均值到槳葉損傷的線性觀測方程為:

其中:yk為槳葉損傷第k次觀測值，H1為診斷矩陣，h2為常數(shù)項。設(shè)得到關(guān)于槳葉損傷及殘差均值的實驗樣本，則系數(shù)H1，h2可確定方法為:

式中:

由于傳感器測量誤差、Volterra預(yù)測模型在不同區(qū)域內(nèi)的精度變化、靈敏度在相軌跡上的起伏誤差等隨機干擾影響，由式(9)確定的損傷觀測中通常仍有較大起伏，槳葉損傷狀態(tài)的估計在技術(shù)實質(zhì)上是一種優(yōu)化濾波的過程?？紤]到損傷演化過程(式(2))的非線性，本文選擇Unscented濾波器對損傷觀測數(shù)據(jù)處理。與傳統(tǒng)濾波器(如擴展Kalman濾波器)相比，該濾波器采用Unscented變換通過2n+1個西格瑪點及狀態(tài)轉(zhuǎn)移點的聯(lián)合分布對狀態(tài)變量進行估計，可較好反映非線性狀態(tài)轉(zhuǎn)移函數(shù)中的高階距成分，降低分析誤差［11］。Unscented濾波器主要含預(yù)測步和修正步，預(yù)測步用于計算損傷變量均值及方差的一步估計值，計算式為:

修正步用于估計觀測變量均值及方差，并用觀測結(jié)果修正損傷變量的均值及方差，計算式為:

式中:Kk為濾波增益;Rk為觀測噪聲方差。

圖10 采用UKF濾波器估計槳葉損傷流程Fig.10 Flowchart to estimate blade damage using Unscented Filter

圖11 槳尖揮舞位移對損傷的估計Fig.11 Tracking results of damage using blade flaps

3 結(jié)論

本文用信號重構(gòu)技術(shù)與Volterra預(yù)測模型提出新的損傷觀測信號生成方法，該方法利用信號的非線性特性在更高維數(shù)的狀態(tài)空間內(nèi)重構(gòu)旋翼系統(tǒng)的動力學(xué)本質(zhì)，并運用現(xiàn)代濾波技術(shù)在飛行狀態(tài)下跟蹤旋翼的損傷狀況，克服了旋翼可測狀態(tài)變量數(shù)目的限制問題。對結(jié)構(gòu)損傷機理和材料失效理論的研究通常在微觀尺度上進行。但此研究不能闡明損傷過程在整個系統(tǒng)宏觀行為中的表現(xiàn)，不能提供系統(tǒng)級損傷跟蹤手段。而源自混沌分析理論的相空間重構(gòu)技術(shù)，可從不同時空角度觀察系統(tǒng)損傷行為特性。研究表明，只要時間尺度及空間維數(shù)選擇合適，即可通過相軌跡光滑移動了解損傷的累積過程。

［1］ Azzam H，Andrew M J.The use of math-dynamic model to aid the development of integrated health and usage monitoring［J］.Journal of Aerospace Engineering，1992，206(1):71-96.

［2］ Ganguli R，Chopra I，Haas D J.Detection of helicopter rotor system simulated faults using neural networks［J］.Journal of American Helicopter Soc，1997，42(2):161-171.

［3］Pawar P M，Ganguli R.Fuzzy logic based health monitoring and residual life prediction for composite rotor［J］.Journal of Aircraft，2007，44(3):981-995.

［4］Pawar P M，Ganguli R.Helicopter rotor health monitoring-a review［J］.Journal of Aerospace Engineering，2007，221(1):631-647.

［5］Johnson W.Comprehensive analytical model for rotorcraft aerodynamics and dynamics volumes ii:components theory［M］.California:Johnson Aeronautics，Palo Alto，2007:280-283.

［6］ Miles T J，Lucas M.Torsional and bending vibration measurements on rotors using laser technology［J］.J.Sound Vibr.，1999，266(3):1441-1467.

［7］Pawar P M，Ganguli R.On the effect of progressive damage on composite helicopter rotor system behavior［J］.Composite Structures，2007，78(3):410-423.

［8］ Robert C H.Chaos and nonlinear dynamics［M］.Oxford University Press，2002.

［9］陳 果.非線性時間序列的動力學(xué)混沌特征自動提取技術(shù)［J］.航空動力學(xué)報，2007，22(1):1-6.CHEN Guo.Dynamic chaos features auto extracting technique of nonlinear time series［J］.Journal of Aerospace Power，2007，22(1):1-6.

［10］張家樹，肖先賜.混沌時間序列的 Volterra自適應(yīng)預(yù)測［J］.物理學(xué)報，2000，49(3):404-408.

ZHANG Jia-shu，XIAO Xian-ci.Predicting low dimensional chaotic time series using volterra adaptive filer［J］.Acta Physic Asinic，2000，49(3):404-408.

［11］ Julier S J，Uhlmann J K.Unscented filtering and nonlinear estimation［J］.Proc.IEEE，2004，92(3):401-422.

［12］ Chelidze D，Liu M.Dynamical systems approach to fatigue damage identification［J］.J Sound and Vibration，2005，281(3):887-904.