基于混沌吸引子的飛輪故障檢測

李磊， 高永明， 吳止鍰

(航天工程大學航天信息學院, 北京 101416)

隨著航天技術的發展，航天器的故障診斷技術越來越受到重視。姿態控制系統是航天器的重要子系統，航天器的功能發揮依賴于特定的姿態，姿態控制系統一旦發生故障，將對航天器造成致命的影響。很多故障在發生之前會表現為早期故障征兆[1]，早期故障的檢測對于故障預報有著重要意義。

針對飛輪的故障診斷主要有2類方法[2]，包括基于解析模型的方法和基于知識的方法。基于解析模型的方法包括觀測器[3-4]和濾波器[5-6]。這一類方法需要建立系統的數學模型，模型精度決定故障診斷效果，因此這種方法限制條件較多，可移植性差。基于知識的方法包括神經網絡[7-8]和支持向量機[9-10]。這一類方法無需建立系統的精確數學模型，通過歷史數據的訓練能夠得到足夠精度的擬合模型，因此該方法具有學習能力，但是其參數調整十分復雜，且模型的可解釋性差。為了避免復雜的數學建模，文獻[11-13]基于數據驅動實現了飛輪的故障診斷，但是也存在計算復雜、數據量要求高等不足。

混沌性是一種廣泛存在于非線性系統中的特性[14]，相空間重構理論[15]表明可以利用較少的變量刻畫復雜系統的特性，因此混沌性可以作為系統特征用于故障診斷。基于混沌的故障診斷在電力、機械領域已經取得了應用[16-17]，在航天領域亦有少量研究[18-19]，具有對微弱故障的敏感性和噪聲的免疫力。

對于飛輪這一類物理模型復雜的非線性系統，本文提出一種數據驅動的故障檢測方法；由于飛輪的混沌性在控制算法的作用下不顯著，構造一個由飛輪參數和輔助函數組成的離散動力系統，該系統可以穩定地產生混沌吸引子；利用混沌吸引子形態在故障情況下的變化檢測故障。

1 飛輪模型

根據文獻[20]給出的ITHACO-A型反作用飛輪的高精度數學模型，建立了基于Simulink的飛輪高精度仿真模型，如圖1所示。

圖1 高精度飛輪仿真模型Fig.1 High-accuracy simulation model of flywheel

圖1中相關參數如下：Vc為等效控制電壓；τz為輸出力矩；Vbus為母線電壓；ω為飛輪轉速；ωd為驅動帶寬；Im為電機電流；Ibus為母線電流；B、C為常數；Hb、Hf和Hz均為Heaviside step函數；參考文獻[13]設置仿真參數，具體取值見表1。

飛輪的常見故障模式有：卡死故障、空轉故障、摩擦增大故障、增益下降故障和跳變故障。按照故障隨時間的變化可分為突變故障和緩變故障2種。緩變故障在早期通常會表現為幅度較小的性能變化，其中摩擦增大故障和增益下降故障是較常見的早期故障。摩擦增大故障主要表現為飛輪的輸出力矩與期望控制力矩差異增大，其故障模型可等效為各類摩擦系數變大；增益下降故障主要表現為飛輪輸出力矩相對期望控制力矩比例減小，其故障模型可等效為電機驅動增益Gd減小。

表1 飛輪仿真模型參數設置

根據圖1建立的飛輪模型可以得到如下的飛輪狀態空間模型：

(1)

式中：Vvc為飛輪輸入的控制電壓；ψ1(Im,ω)、ψ2(ω)和ψ3(ω)分別為電動勢力矩限制、軸承摩擦力與力矩干擾和電機轉速限制這3個非線性過程。式(1)反映了電機電流是表征系統的重要變量，系統狀態的改變必然伴隨著電機電流的變化，因此本文選用電機電流Im作為分析的變量。

2 故障檢測方法

基于混沌的故障檢測方法通常利用參數時間序列自身的混沌特征進行檢測，如最大Lyapunov指數、分形維數、盒維數等。該方法的前提條件是待檢測系統處于混沌狀態。飛輪系統由于處在控制算法作用之下，其混沌特性被抑制，并不能總是處在混沌狀態中，因此通常的基于混沌的故障檢測方法不適用于飛輪系統的故障檢測。本文采取構造混沌離散動力系統的方法進行飛輪系統的故障檢測。

2.1 離散動力系統的構造方法

以第1節所述模型進行混沌性分析，當飛輪保持在任意穩定工況下時，計算Im不同時間數據序列的最大Lyapunov指數λ，結果見圖2，n為采樣點數。其值并不總是正的，即飛輪系統的混沌性不是穩定出現的，也意味著不能直接使用系統的混沌特性。因此需要構造一個存在穩定的混沌特性的系統。

能夠產生混沌的動力系統模型有很多，常見的如Lorenz系統、Logistic系統和Henon系統等。但是將飛輪系統的數學模型變換為類似的混沌系統形式是比較困難的，因此考慮其他形式的構造方法。文獻[21]的研究表明對于圖像函數而言，存在滿足一定條件的輔助函數使得構成的離散動力系統出現混沌吸引子。基于這一結論構造如下離散動力系統：

(2)

式中：k為正弦曲面函數參數；相空間重構灰度圖像g(x,y)由待分析變量Im經過相空間重構生成；f(x,y)和g(x,y)的定義域與值域應調整一致。選擇該輔助函數的原因是在參數k的控制下能夠形成不同振蕩程度、不同梯度以及充滿整個定義域的曲面，這樣的輔助曲面更容易產生混沌。

圖2 Im的最大Lyapunov指數Fig.2 The largest Lyapunov exponent of Im

2.2 離散動力系統參數的選擇

構造g(x,y)有2個目的：一是使其符合灰度圖像函數的特點；二是能反映原數據序列的特征。相空間重構理論證明了存在一個合適的嵌入維m，如果滿足m>2d+1，d為動力系統的維數，那么在這個m維的相空間中可以恢復出原系統的吸引子。通過選擇不同的嵌入維m和時延τ組合，可以生成不同大小的相空間矩陣，對應不同尺寸的灰度圖像函數。因此采用相空間重構法生成所需函數。

對于一個單變量時間序列{xi}，如果嵌入維和時延為m和τ，則可按如式(3)形式構造相空間[15]:

Xi=[xi,xi+τ,L,xi+(m-1)τ]

(3)

式中：i=1，2，…，L；L=N-(m-1)τ，N為序列長度。

對于一般的系統而言，不失一般性的，時延τ通常可以選擇為1或者2，具有簡單性、比較強的可操作性和實用性。嵌入維m的選擇需要考慮灰度圖像的特點，本文選擇2的冪數256、128或64等。不同維數對結果的影響會在3.3節予以討論。

為確保構造離散動力系統能夠大概率的產生混沌吸引子，需要選擇合適的輔助函數f(x,y)參數k。為此給出了在m=256的條件下，參數k變化下的分岔圖和系統最大Lyapunov指數圖，分別如圖3和圖4所示。

圖3 離散動力系統的分岔圖Fig.3 Bifurcation diagram of discrete dynamic system

圖4 離散動力系統的最大Lyapunov指數Fig.4 The largest Lyapunov exponent of discrete dynamic system

如圖3所示，當參數k大于3時，系統進入混沌狀態。對比圖4，在同樣的k值下絕大部分最大Lyapunov指數為正，證明了在這一k值區間中系統會產生混沌吸引子。

2.3 基于混沌吸引子的故障檢測原理分析

混沌系統具有以下性質：對初值敏感性、有界性、遍歷性、內隨機性和分維性等。主要利用前3種性質進行故障檢測。混沌系統對初值敏感性表明，對于一個混沌系統，不同的初值會產生不同的響應；有界性表明系統的運動軌跡始終局限在一個確定的區域內；遍歷性表明系統運動軌跡在吸引域內是各態遍歷的。

由于混沌吸引子的有界性和遍歷性，使用2.1與2.2節構造的離散動力系統生成近似吸引子時只要迭代次數足夠，則無論選擇哪個初始點，最終的吸引子形狀是一致的。圖5(a)為初始點(65，90)開始迭代300次后的吸引子，圖5(b)為初始點(90，65)開始迭代300次后的吸引子，二者形狀基本一致。

當飛輪系統發生故障時，參數變化規律發生改變，構造的離散動力系統特性也隨之改變，于是迭代生成的混沌吸引子形狀發生變化。圖6為30%幅度的飛輪電機增益變小故障下吸引子圖像，與圖5(a)和圖5(b)相比吸引子的形狀變化顯著。利用這一變化可以進行故障檢測。

圖5 初始點(65,90)和(90，65)的混沌吸引子Fig.5 Chaotic attractor with initial points(65,90) and (90,65)

為利用混沌吸引子的變化進行故障檢測，需要提取描述混沌吸引子的特征。本文構造產生的混沌吸引子實際是一個二位平面函數，因此可以采用rodon變換將其映射為一維數組，通過對變換后的一維數組進行多項式擬合，即可得到混沌吸引子的特征多項式。正常情況下的混沌吸引子特征應是近似相同的，而故障情況下的混沌吸引子特征與正常情況應有較大區別，因此計算不同混沌吸引子特征多項式的相關系數進行故障檢測。

混沌系統對初值的敏感性使得基于混沌吸引子的故障檢測方法對故障十分敏感，這一點會在3.2節的仿真中進行驗證。

圖6 初始點(65,90)，Kgd=0.7時的混沌吸引子Fig.6 Chaotic attractor with initial point(65,90) and Kgd=0.7

2.4 基于混沌吸引子的故障檢測方法流程

2.3節的分析中論述了混沌吸引子與故障間的關系，以此為依據設計如下的故障檢測方法：

1) 數據預處理。按所需的灰度圖像尺寸采集相應長度的數據序列，對采集到的數據進行歸一化處理使數據分布在[0,1]區間內，將數據值轉換為灰度值。

2) 生成混沌吸引子。隨機選擇初始點，以構造的離散動力系統進行一定次數的迭代計算生成混沌吸引子。

3) 計算混沌吸引子特征。對生成的混沌吸引子矩陣進行radon變換，對變換后的數據進行擬合，以擬合的多項式生成固定長度的特征向量。

4) 故障檢測。選取平穩工況下的混沌吸引子作為標準，計算待檢測序列的混沌吸引子與其相關系數，相關系數低于閾值則認為檢測到故障的發生。

3 數值仿真分析

3.1 閉環飛輪仿真系統

在文獻[13]的仿真條件基礎上，增加了不同工況下的故障檢測仿真。搭建的閉環飛輪仿真系統如圖 7所示。

(4)

考慮實際的姿態控制系統僅少部分時間處在大角度機動狀態下，大部分時間為穩定姿態，修正空間干擾力矩以及軌道轉速的影響，因此本文指定的飛輪工作狀態為抵抗空間干擾力矩作用下的姿態保持過程，設置目標轉速為常值疊加干擾力矩對應的速度變化：

(5)

式中：ωc=4.5 rad/s；ω0=0.02 rad/s；A=0.14。

設置3種測量噪聲條件，信噪比分別為75 dB(標稱)、55 dB(中等強度)和45 dB(高強度),其中“標稱”噪聲是指目前器件的典型噪聲水平。設置了電機增益變小和摩擦力矩增大2種故障。每種噪聲條件下2種故障各有6組數據，其中1組為正常數據；5組為故障數據。設置電機增益變小故障系數Kgd分別為0.7、0.82、0.91、0.94和0.97，即故障參數的偏差幅度分別為正常值的30%、18%、9%、6%和3%；摩擦力矩故障系數Kτ分別為1.3、1.21、1.09、1.06和1.03，即故障參數的偏差幅度分別為正常值的30%、18%、9%、6%和3%。

圖7 閉環飛輪控制系統框圖Fig.7 Block diagram of flywheel with closed-loop control system

故障檢測的閾值根據正常數據產生的混沌吸引子相似度選取，可以選擇最大值或是平均值，本文選取最大值為閾值。

系統仿真長度為400 s，步長為0.01 s，設置故障時刻為200 s，單次仿真共40 000 點數據。仿真過程中電流和飛輪轉速初值均為零。采用的分析量為電機電流Im，由于系統初始時刻有一個調節過程，因此本文只從50 s后的穩態開始進行分析，即每次仿真采用50～400 s間的數據，共35 000個采樣點。

3.2 飛輪系統的故障檢測

3.2.1 電機增益變小故障

標稱噪聲(75 dB)條件下的5段故障數據均為50～400 s間的采樣點，每段為35 000個采樣點，如圖 8所示。微小幅度下的故障與正常數據的偏差很小，基本掩蓋在了噪聲范圍內，當故障幅度逐漸增大時，故障數據與正常數據有了明顯的偏差，但是經過一段時間后，在反饋控制的作用下系統重新恢復穩定。

按2.4節所述方法對采集的數據進行計算，其中輔助函數參數k取值為3.141 6，嵌入維m取值為256。

標稱噪聲(75 dB)條件下不同故障幅度的檢測結果如圖9所示，kre為相差系數。隨著故障幅值的增大，其與正常數據的相似度越來越低，以0.8為閾值的情況下，對3%幅度(即Kgd=0.97)的故障依然響應靈敏，該方法對電機增益變大故障很敏感。同時，由于控制目標的波動，圖9也表明了飛輪系統的平穩狀態特征與電機增益故障特征是與控制目標的幅值無關，因此該方法可以應用在不同控制目標的條件下。

圖8 飛輪電機增益變小故障的電機電流Fig.8 Current of flywheel under fault of motor gain decrease

考察本文方法對噪聲的魯棒性，以3種噪聲條件下故障系數Kgd=0.7的故障檢測為例，檢測結果如圖10所示。相同的故障在不同噪聲條件下的檢測效果是一致的，證明本文方法對噪聲具有魯棒性。

圖9 75 dB噪聲下電機增益變小故障檢測結果Fig.9 Fault detection result of motor gain decrease under 75 dB noise

圖10 不同噪聲條件下電機增益變小故障檢測結果(Kgd=0.7)Fig.10 Fault detection result of motor gain decrease under different noises (Kgd=0.7)

3.2.2 摩擦力矩增大故障

標稱噪聲(75 dB)條件下的5段故障數據均為50～400 s間的采樣點，每段為35 000個采樣點。如圖11所示，微小幅度下的故障與正常數據的偏差很小，大部分處在噪聲范圍內，當故障幅度逐漸增大時，故障數據與正常數據有了明顯的偏差，由于摩擦力增大故障是通過轉速間接影響電機電流的，因此在經過一段時間后，在反饋控制的作用下系統重新恢復穩定時的電流值會變大。

按2.4節所述方法對采集的數據進行計算，其中輔助函數參數k取值為3.141 6，嵌入維m取值為256。

標稱噪聲(75 dB)條件下不同故障幅度的檢測結果如圖12所示。在閾值為0.8時能夠檢出9%幅度(即Kτ=1.09)的故障，但對于更小幅度故障的檢測區分度不如電機增益變大故障的檢測明顯。同時，由于控制目標的波動，圖12也表明了飛輪系統的平穩狀態特征與電機增益故障特征是與控制目標的幅值無關，因此該方法可以應用在不同控制目標的條件下。

考察本文方法對噪聲的魯棒性，以3種噪聲條件下故障系數Kτ=1.3的故障檢測為例，檢測結果如圖13所示。相同的故障在不同噪聲條件下的檢測效果是一致的，證明本文方法對噪聲具有魯棒性。但是與電機增益故障相比，摩擦力矩故障對噪聲更敏感。

圖11 摩擦力矩增大故障的的電機電流Fig.11 Current of motor in fault of motor friction moment increase

圖12 75 dB噪聲下摩擦力矩增大故障檢測結果Fig.12 Fault detection result of motor friction moment increase under 75 dB noise

3.3 嵌入維對故障檢測能力的影響

以上仿真選擇的嵌入維均為256，這是因為通常灰度圖像都具有256個級別的灰度值，在256×256的尺寸上無需調整灰度值，方便構造混沌吸引子。但是這樣對數據長度的要求相對較高，需要766長度的數據。下面討論不同嵌入維對檢測效果的影響。

為方便計算，選擇的嵌入維分別為256、128、64和32；所需的數據長度分別為766、382、190和94。以標稱噪聲(75 dB)條件下電機增益故障和摩擦力矩故障為例進行分析，電機增益故障系數Kgd為0.91和0.7；摩擦力矩故障系數Kτ為1.09和1.3。結果如圖14所示。

圖13 不同噪聲條件下摩擦力矩增大故障檢測結果(Kτ=1.3)Fig.13 Fault detection result of motor friction moment increase under different noises (Kτ=1.3)

圖14(a)、(c)的結果顯示，隨著嵌入維數的減小，該檢測方法對故障的敏感度上升，當嵌入維為64時，對于幅度為9%的微小故障識別能力顯著，當嵌入維為32時噪聲的干擾掩蓋了故障。圖14(b)、(d)的結果則顯示，對與故障幅度較大的故障，合理的嵌入維對檢測結果影響不大，但是檢測窗口會隨著嵌入維的減小而減小。

圖14 不同m下電機摩擦力矩增大故障檢測結果Fig.14 Fault detection result of motor friction moment increase with different m

圖15 不同m下電機故障電流均值Fig.15 Mean of fault motor current under different m

4 結 論

1) 仿真結果表明基于混沌吸引子的故障檢測方法能夠很好地檢測出3%幅度的電機增益故障和9%幅度的摩擦力矩故障，由于緩變故障在早期通常會表現為幅度較小的性能變化，因此該方法可以用于飛輪的早期故障檢測。

2) 飛輪系統的數據使用曲面迭代的方法可以產生混沌吸引子，且該吸引子與迭代初始點無關；故障情況下的混沌吸引子形狀會產生明顯的改變。

3) 基于混沌吸引子的故障檢測方法無需對系統進行精確建模，只需要健康數據即可進行故障檢測，同時計算過程大部分為迭代計算，計算量小，適合用于實時故障檢測。

4) 嵌入維數會影響方法對故障的敏感度和檢測窗口。隨著嵌入維數的減小，方法需要的數據減少，對故障的敏感程度增加，但易受噪聲干擾。

5) 基于混沌吸引子的故障檢測方法是一種基于數據的方法，可移植性強，可以推廣應用到航天器其他系統、部件的故障檢測。