李磊， 高永明， 吳止鍰

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

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

針對飛輪的故障診斷主要有2類方法[2]，包括基于解析模型的方法和基于知識的方法?；诮馕瞿Ｐ偷姆椒òㄓ^測器[3-4]和濾波器[5-6]。這一類方法需要建立系統(tǒng)的數(shù)學(xué)模型，模型精度決定故障診斷效果，因此這種方法限制條件較多，可移植性差?；谥R的方法包括神經(jīng)網(wǎng)絡(luò)[7-8]和支持向量機[9-10]。這一類方法無需建立系統(tǒng)的精確數(shù)學(xué)模型，通過歷史數(shù)據(jù)的訓(xùn)練能夠得到足夠精度的擬合模型，因此該方法具有學(xué)習(xí)能力，但是其參數(shù)調(diào)整十分復(fù)雜，且模型的可解釋性差。為了避免復(fù)雜的數(shù)學(xué)建模，文獻[11-13]基于數(shù)據(jù)驅(qū)動實現(xiàn)了飛輪的故障診斷，但是也存在計算復(fù)雜、數(shù)據(jù)量要求高等不足。

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

對于飛輪這一類物理模型復(fù)雜的非線性系統(tǒng)，本文提出一種數(shù)據(jù)驅(qū)動的故障檢測方法；由于飛輪的混沌性在控制算法的作用下不顯著，構(gòu)造一個由飛輪參數(shù)和輔助函數(shù)組成的離散動力系統(tǒng)，該系統(tǒng)可以穩(wěn)定地產(chǎn)生混沌吸引子；利用混沌吸引子形態(tài)在故障情況下的變化檢測故障。

1 飛輪模型

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

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

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

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

表1 飛輪仿真模型參數(shù)設(shè)置

根據(jù)圖1建立的飛輪模型可以得到如下的飛輪狀態(tài)空間模型：

(1)

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

2 故障檢測方法

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

2.1 離散動力系統(tǒng)的構(gòu)造方法

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

能夠產(chǎn)生混沌的動力系統(tǒng)模型有很多，常見的如Lorenz系統(tǒng)、Logistic系統(tǒng)和Henon系統(tǒng)等。但是將飛輪系統(tǒng)的數(shù)學(xué)模型變換為類似的混沌系統(tǒng)形式是比較困難的，因此考慮其他形式的構(gòu)造方法。文獻[21]的研究表明對于圖像函數(shù)而言，存在滿足一定條件的輔助函數(shù)使得構(gòu)成的離散動力系統(tǒng)出現(xiàn)混沌吸引子?；谶@一結(jié)論構(gòu)造如下離散動力系統(tǒng)：

(2)

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

圖2 Im的最大Lyapunov指數(shù)Fig.2 The largest Lyapunov exponent of Im

2.2 離散動力系統(tǒng)參數(shù)的選擇

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

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

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

(3)

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

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

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

圖3 離散動力系統(tǒng)的分岔圖Fig.3 Bifurcation diagram of discrete dynamic system

圖4 離散動力系統(tǒng)的最大Lyapunov指數(shù)Fig.4 The largest Lyapunov exponent of discrete dynamic system

如圖3所示，當(dāng)參數(shù)k大于3時，系統(tǒng)進入混沌狀態(tài)。對比圖4，在同樣的k值下絕大部分最大Lyapunov指數(shù)為正，證明了在這一k值區(qū)間中系統(tǒng)會產(chǎn)生混沌吸引子。

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

混沌系統(tǒng)具有以下性質(zhì)：對初值敏感性、有界性、遍歷性、內(nèi)隨機性和分維性等。主要利用前3種性質(zhì)進行故障檢測?；煦缦到y(tǒng)對初值敏感性表明，對于一個混沌系統(tǒng)，不同的初值會產(chǎn)生不同的響應(yīng)；有界性表明系統(tǒng)的運動軌跡始終局限在一個確定的區(qū)域內(nèi)；遍歷性表明系統(tǒng)運動軌跡在吸引域內(nèi)是各態(tài)遍歷的。

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

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

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

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

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

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

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

2.3節(jié)的分析中論述了混沌吸引子與故障間的關(guān)系，以此為依據(jù)設(shè)計如下的故障檢測方法：

1) 數(shù)據(jù)預(yù)處理。按所需的灰度圖像尺寸采集相應(yīng)長度的數(shù)據(jù)序列，對采集到的數(shù)據(jù)進行歸一化處理使數(shù)據(jù)分布在[0,1]區(qū)間內(nèi)，將數(shù)據(jù)值轉(zhuǎn)換為灰度值。

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

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

4) 故障檢測。選取平穩(wěn)工況下的混沌吸引子作為標(biāo)準(zhǔn)，計算待檢測序列的混沌吸引子與其相關(guān)系數(shù)，相關(guān)系數(shù)低于閾值則認為檢測到故障的發(fā)生。

3 數(shù)值仿真分析

3.1 閉環(huán)飛輪仿真系統(tǒng)

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

(4)

考慮實際的姿態(tài)控制系統(tǒng)僅少部分時間處在大角度機動狀態(tài)下，大部分時間為穩(wěn)定姿態(tài)，修正空間干擾力矩以及軌道轉(zhuǎn)速的影響，因此本文指定的飛輪工作狀態(tài)為抵抗空間干擾力矩作用下的姿態(tài)保持過程，設(shè)置目標(biāo)轉(zhuǎn)速為常值疊加干擾力矩對應(yīng)的速度變化：

(5)

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

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

圖7 閉環(huán)飛輪控制系統(tǒng)框圖Fig.7 Block diagram of flywheel with closed-loop control system

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

系統(tǒng)仿真長度為400 s，步長為0.01 s，設(shè)置故障時刻為200 s，單次仿真共40 000 點數(shù)據(jù)。仿真過程中電流和飛輪轉(zhuǎn)速初值均為零。采用的分析量為電機電流Im，由于系統(tǒng)初始時刻有一個調(diào)節(jié)過程，因此本文只從50 s后的穩(wěn)態(tài)開始進行分析，即每次仿真采用50～400 s間的數(shù)據(jù)，共35 000個采樣點。

3.2 飛輪系統(tǒng)的故障檢測

3.2.1 電機增益變小故障

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

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

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

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

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

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

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

3.2.2 摩擦力矩增大故障

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

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

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

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

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

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

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

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

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

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

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

圖14 不同m下電機摩擦力矩增大故障檢測結(jié)果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 結(jié) 論

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

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

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

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

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