基于故障分辨能力的火箭發(fā)動機測量參數(shù)選擇方法

張效溥，任楓，徐鵬里，李志敏，宿彩虹

上海宇航系統(tǒng)工程研究所，上海 201109

液體火箭發(fā)動機是火箭動力系統(tǒng)的關(guān)鍵組成部分，其工作狀態(tài)直接關(guān)系到發(fā)射任務(wù)的成敗。目前，隨著載人登月、深空探測等攻關(guān)項目的開展，為滿足“一重故障正常工作，兩重故障保證安全”的功能要求，需協(xié)同開發(fā)故障監(jiān)測、健康管理、飛行重構(gòu)等功能，這對測控系統(tǒng)對故障的分辨能力和傳感器的可靠性提出了巨大的挑戰(zhàn)。例如，對某現(xiàn)役飛行器發(fā)動機開展故障分析時，發(fā)現(xiàn)故障特征可能由5 種不同的失效形式造成，無法定位到具體單機產(chǎn)品，嚴重影響故障排查進展；在某現(xiàn)役型號火箭飛行重構(gòu)策略的制定中，也出現(xiàn)了在當(dāng)前測量方案下末子級發(fā)動機大量故障特征相似程度極高、難以進行診斷的問題。依據(jù)研制經(jīng)驗，傳感器布置的數(shù)量、安裝位置、安裝形式對遙測數(shù)據(jù)的獲取和發(fā)動機的可靠性都有著直接的影響。若傳感器布置太多，可能造成測量特征的冗余和浪費，加大了總裝總測的難度，且降低了發(fā)動機的可靠性。以往飛行試驗就曾出現(xiàn)因在高壓容腔安裝測壓管接頭造成的角焊縫撕裂的故障模式；相反，盲目地減少傳感器數(shù)量，會使大量故障的特征趨于一致，造成責(zé)任界面的不明確、“舉一反三”工作難開展。如何選取發(fā)動機的測量參數(shù)并布置傳感器，使發(fā)動機保證較高可靠性的同時分辨更多種類的性能故障，是一個關(guān)系到發(fā)動機設(shè)計和使用的關(guān)鍵問題。

測量特征的選擇是一個典型的最優(yōu)化問題，通過不斷調(diào)整特征子集或者測點的空間屬性改善關(guān)心的測量性能。Maul 等［1］以飛機系統(tǒng)為例，提出了傳感器布置的4 類評價指標：性能的可觀測性、故障的檢測能力、可靠性、費用。目前國內(nèi)外相關(guān)研究大多圍繞這4 類指標展開。鄭帥［2］、Jung［3］等對以不可測油量和姿態(tài)誤差為優(yōu)化目標，對飛機油箱內(nèi)的液位傳感器的安裝位置和角度開展了優(yōu)化；Pereira 等［4］將7 種性能指標作為優(yōu)化目標，對直升機葉片測量方案開展優(yōu)化；Omata 等［5］采用主成分分析對火箭發(fā)動機不同位置泄漏故障進行降維，并采用貪婪策略搜索更容易檢測出泄漏故障的傳感器排布；徐敏強等［6］采用符號有向圖完成火箭發(fā)動機故障專家知識的推斷，并提出一種提高故障檢測覆蓋率的方法；楊冬健等［7］梳理了航空液壓系統(tǒng)故障因果矩陣，并提出了一種提高故障隔離率、降低虛警率的傳感器選擇方法；張笑華［8］、Li［9］、Liu［10］等將傳感器優(yōu)化技術(shù)應(yīng)用于結(jié)構(gòu)健康監(jiān)測，通過優(yōu)化加速度傳感器空間位置，加強對結(jié)構(gòu)損傷的測量靈敏度；Kong 等［11］基于專家知識，采用離散粒子群算法對液壓控制系統(tǒng)的傳感器排布開展優(yōu)化，提高故障檢測的實時性；Li 和Der Kiureghian［12］采用最大期望效用理論和貝葉斯線性模型來進行魯棒傳感器的布置，以實現(xiàn)運行模態(tài)的辨識；Li等［13］將航空發(fā)動機進氣道傳感器優(yōu)化任務(wù)轉(zhuǎn)化為同時含1、2 范數(shù)懲罰項的最優(yōu)化問題，并提出了一種收斂性較好的求解方法。此外，部分研究人員嘗試將各領(lǐng)域的傳感器優(yōu)化問題形成規(guī)范化的優(yōu)化算法框架。例如，B?achowski 等［14］提出一種高效的動載荷下大型建筑傳感器布置優(yōu)化框架；Taravatrooy 等［15］在考慮了不確定性的前提下，提出了改善供水管網(wǎng)系統(tǒng)泄漏檢測能力的壓力傳感器優(yōu)化布置算法框架；Morlier 等［16］基于Kriging 方法建立傳感器位置信息與MAC 系數(shù)間的代理模型，形成傳感器布置高效全局優(yōu)化方案。

從上述調(diào)研可以看出，當(dāng)前應(yīng)用于液壓、氣動系統(tǒng)傳感器優(yōu)化布置主要基于專家知識或者信息論，但是這2 類方法對火箭發(fā)動機傳感器優(yōu)化均存在一定的局限性?；趯＜抑R的方法更多地關(guān)注故障的“覆蓋性”，即“是否能測到”而不是“是否能分清”，沒有形成量化的故障特征庫，對故障間的距離度量考慮不足，難以解決“不同類型故障相似度較高”這一困擾；同時，基于專家知識的傳感器優(yōu)化方法非常依賴于領(lǐng)域知識的正確性，而火箭發(fā)動機作為復(fù)雜閉環(huán)系統(tǒng)，故障種類眾多，難以根據(jù)工程經(jīng)驗開展故障分析；基于信息論的傳感器優(yōu)化方法一般采用相關(guān)系數(shù)、冗余指標、互信息系數(shù)等作為優(yōu)化目標。這些指標在一定程度上可以反映出傳感器布置方案的優(yōu)劣，但是不能直接表征可區(qū)分故障的數(shù)量和類別。同時，現(xiàn)有方法鮮有考慮制造公差、測量噪聲下的魯棒性和傳感器安裝對系統(tǒng)的可靠性。

基于上述現(xiàn)狀，提出了一種基于模型的、面向液體火箭發(fā)動機的測點組合對故障分類能力的評估方法以及優(yōu)化方法。為克服火箭發(fā)動機故障專家知識提取較為困難的問題，基于發(fā)動機靜態(tài)性能仿真模型構(gòu)建發(fā)動機故障特征庫；采用凝聚層次聚類算法評估任一測點排布能識別的故障種類，并結(jié)合測點組合的抗噪性能指標、測點組合的風(fēng)險量化指標共同構(gòu)成測點綜合評價指標；采用二進制多目標粒子群算法對測點組合進行優(yōu)化，獲得測點組合的Pareto 解集，并挑選優(yōu)化后的測點組合。該方法可以較好地平衡系統(tǒng)可靠性、故障的可分性和故障識別的魯棒性，對其他復(fù)雜、閉環(huán)的動力系統(tǒng)的測點選取亦有良好的應(yīng)用前景。

1 發(fā)動機靜態(tài)特性建模及故障庫的建立

1.1 液體火箭發(fā)動機數(shù)學(xué)模型

為建立液體火箭發(fā)動機的靜態(tài)模型，應(yīng)先確定發(fā)動機各組件的靜態(tài)特性方程，再代入系統(tǒng)約束方程進行求解。從發(fā)動機結(jié)構(gòu)來看，發(fā)動機由管路、閥門、燃燒室、渦輪等組件組成；從發(fā)動機系統(tǒng)設(shè)計來看，主要由主系統(tǒng)和副系統(tǒng)構(gòu)成，且滿足穩(wěn)態(tài)壓力、流量和功率3 大方程。下面介紹各關(guān)鍵單機數(shù)學(xué)模型［17］。

1.1.1 渦輪泵模型

渦輪泵是泵壓式發(fā)動機的核心組件之一，其作用是將燃氣的焓轉(zhuǎn)化為機械功，并對推進劑進行增壓，使其流入推力室。渦輪泵由1 個軸流式氣動渦輪和2 個離心式推進劑泵組成。其中，渦輪功率、渦輪效率分別為

式中：Wt為渦輪功率；ηt為渦輪效率；γ為氣體絕熱系數(shù)；Rg為燃氣的氣體常數(shù)；T*為燃氣總溫；pin、pout分別為渦輪上、下游壓力；qf為燃氣流量；n為轉(zhuǎn)速；at、bt、ct均為渦輪效率常數(shù)。

離心泵的功率、揚程與效率特性分別為

式中：Wp、Δpb、ηb分別為泵的功率、揚程和效率；qv為推進劑體積流量；n0為額定轉(zhuǎn)速；ah、bh、ch均為泵揚程的渦輪泵特性常數(shù)；ap、bp、cp均為泵效率常數(shù)。

1.1.2 燃燒組件模型

燃燒組件的作用是將雙組元推進劑的化學(xué)能轉(zhuǎn)化為燃氣的內(nèi)能，形成較高的室壓、室溫，再經(jīng)過噴嘴、噴管元件形成流量。燃燒組件的室溫、室壓和流量的關(guān)系為

式中：c*為燃氣特征速度；pc、qc、At、ηc分別為燃燒組件的室壓、流量、喉徑和燃燒效率；燃氣組分和溫度一般采用等焓約束下的最小自由能法求解，其數(shù)學(xué)模型為

其中：ni、hi和si為第i種燃燒產(chǎn)物的分子量、比焓和比熵；Li，j、Ltot，j分別為第i種燃氣中第j種元素原子的數(shù)量和第j種元素原子總數(shù)；T為燃氣溫度；Mi、mtot分別為第i種組分的分子量和燃燒室內(nèi)物質(zhì)總質(zhì)量；Hin、Hout分別為燃燒前后物質(zhì)的總焓。該問題是含等式約束的最優(yōu)化問題，一般采用Lagrange 乘子法求解。本文為減少計算時間，事先在一定混合比、室壓設(shè)計范圍內(nèi)采樣和計算，采用響應(yīng)面法求解燃氣組分和溫度對混合比、室壓的代理模型，再代入系統(tǒng)模型求解。

對于后端連接拉瓦爾噴管的燃燒組件，需補充推力和比沖的數(shù)學(xué)方程，其表達式為

式中：CFv是噴管的推力系數(shù)；I為噴管的比沖；F為噴管產(chǎn)生的推力。

1.1.3 節(jié)流元件模型

液體火箭發(fā)動機節(jié)流元件有節(jié)流圈、音速噴嘴和汽蝕管3 種。其數(shù)學(xué)特性和功能存在一定的差異。節(jié)流圈起到調(diào)節(jié)推進劑流量的作用，其流量特性為

式中：Qlz為液路節(jié)流圈的流量；Cd為節(jié)流系數(shù)；Ao為節(jié)流孔面積；ρ為推進劑密度。汽蝕管與普通節(jié)流圈相比，具有隔絕下游壓力波動的功能。當(dāng)壓比低于發(fā)生氣蝕的臨界壓比時，汽蝕管退化為節(jié)流圈。其流量特性為

式中：Qqs為氣蝕管的流量；psv為推進劑的飽和蒸氣壓；πcri為發(fā)生汽蝕的臨界壓比。

燃氣噴嘴多用于發(fā)生器后，調(diào)節(jié)渦輪燃氣流量。其流量特性與可壓流體節(jié)流孔類似，可表示為

式中：Qps為音速噴嘴流量；Z為氣體壓縮因子；?為噴嘴的流量系數(shù)。

1.1.4 發(fā)動機系統(tǒng)靜態(tài)模型

按照發(fā)動機結(jié)構(gòu)，形成發(fā)動機熱力學(xué)參數(shù)傳遞的拓撲圖，構(gòu)建系統(tǒng)級平衡方程，即：任意1 個節(jié)點處流量之和為0，任意1 條回路上的壓力變化之和為0，任意1 組同軸的轉(zhuǎn)子功率之和為0。將系統(tǒng)平衡方程聯(lián)立上述單機靜態(tài)模型進行求解即可得到發(fā)動機靜態(tài)參數(shù)。本文采用變步長的Newton-Raphson 方法求解非線性方程組。

1.2 基于模型的故障特征庫的構(gòu)建

發(fā)動機的故障模式按故障類型分類有泄漏、堵塞、工作壓力異常、轉(zhuǎn)子效率下降等；按故障發(fā)生的位置分為主管路故障、副管路故障、渦輪泵故障、閥門故障等。各類故障的定義如表（1）所示。

為全面評估不同故障在各嚴重程度下的發(fā)動機性能變化，每種故障系數(shù)均在相應(yīng)的取值范圍內(nèi)均勻取值30 個，共形成L種故障模式，并計算每種故障模式在這些故障系數(shù)下的故障特征向量。故障特征向量定義為注入故障前后，發(fā)動機m種穩(wěn)態(tài)性能的變化率經(jīng)過規(guī)范化后的向量，其表達式為

式中：yl、yl'為加入第l種故障前、后的發(fā)動機性能參數(shù)；xl為第l種故障造成的性能參數(shù)變化率；xl，norm為歸一化后的性能參數(shù)變化率（1≤l≤L）。將所有故障代入仿真模型進行計算，得到行數(shù)為L、列數(shù)為m的故障特征矩陣Ωall，其定義為

式中：xi，j，norm矩陣第i行、第j列的元素代表第i類故障發(fā)生時第j種發(fā)動機性能的變化率。實際發(fā)動機產(chǎn)品不可能全部測量到所有性能參數(shù)。本文用1 個m維的向量s=［s1，s2，…，sm］來表示發(fā)動機測量參數(shù)的選取情況，其中si取值為0 或1。當(dāng)si為1 時，代表第i種發(fā)動機性能被選為測量參數(shù)；當(dāng)si為0 時，代表第i種特征未被選為測量參數(shù)。則任一測量方案下的故障特征矩陣Ωs為

式中：λs為發(fā)動機測點選擇情況為s時對應(yīng)的故障篩選矩陣。不難看出，該矩陣是Ωall抽取若干列向量組成的新矩陣。

2 考慮故障可分性的測點組合優(yōu)化方法

根據(jù)宇航任務(wù)的特點，發(fā)動機測量參數(shù)選擇應(yīng)滿足以下條件：①盡可能區(qū)分更多故障；② 傳感器便于安裝，且對發(fā)動機系統(tǒng)引入的風(fēng)險盡可能小，可靠性盡可能高；③根據(jù)以往飛行結(jié)果，發(fā)動機產(chǎn)品性能與額定值相比存在一定的隨機誤差，因此測量數(shù)據(jù)應(yīng)具有較好的魯棒性?；诒?，對任意排布方案s，選取凝聚層次聚類的最優(yōu)化劃分數(shù)K作為故障可分性的指標，選取噪聲偏移率e作為魯棒性指標，并基于FMEA 計算測點組合的風(fēng)險指數(shù)R作為可靠性指標，最終形成性能向量［-K（s），e（s），R（s）］作為適應(yīng)度函數(shù)，采用多目標粒子群優(yōu)化算法對測點組合進行尋優(yōu)。

表1 故障模式及定義方式Table 1 Mathematical definition of faults

2.1 故障可分性指標的計算

聚類算法可以較好地表征向量的空間聚集程度，本文基于凝聚層次聚類算法計算故障可分性指標。凝聚的層次聚類是一種自下而上的聚類策略，即把每個單獨的樣本都看作1 個聚類簇，通過不斷合并距離最低、相似度最高的樣本來構(gòu)成新的簇，然后將合并得到的簇看作1 個新的樣本，再次進行合并。如此反復(fù)迭代，直至所有樣本合并為1 類［18］。算法流程為

步驟1設(shè)某測點組合s測量的故障特征矩陣為ΩS。首先將每個行向量都看作1 個單獨的簇Ci，共L組，即Ci=｛xi｝。

步驟2計算每2 個簇之間的距離，形成Lp×Lp的距離矩陣D。其中，dij代表第i個分類簇和第j個分類簇的距離；Lp為第p次合并后的總簇數(shù)。本文選用平均距離法［19］定義簇與簇之間的距離，其定義為

式中：xa、xb為不同簇內(nèi)的向量。

步驟3將距離最近的2 個簇合并成1 個簇，并刪除已合并的2 個簇，如式（21）所示。合并完成后分類簇總數(shù)量減1。

式中：Cnew為合并的簇；Cp、Cq為每步中距離最近的2 個簇。

步驟4重復(fù)步驟2～步驟3，直到所有樣本合并成1 個簇。

步驟5根據(jù)聚類效果度量指標確定最優(yōu)劃分數(shù)K。聚類效果度量指標主要基于簇的緊湊度與分離度，即簇內(nèi)樣本越緊湊、不同簇之間相似度越低，聚類效果越好，越能真實地體現(xiàn)數(shù)據(jù)的聚集情況。常見評價指標的有輪廓系數(shù)、Davies-Bouldin 指數(shù)、Dunn 指數(shù)等［20-21］。本文選取Davies-Bouldin 指數(shù)（DBI）作為度量指標，其定義如式（22）～式（23）所示。

式中：Ip為第p次合并后的DBI。代表第i個和第j個分類簇內(nèi)全部樣本兩兩之間的平均距離，其表達式為

式中：dsingle，i，j為簇內(nèi)第i、第j個向量之間的距離；LCi為簇Ci中的向量數(shù)量。

圖1 為選取某發(fā)動機現(xiàn)有測點對故障表進行聚類得到的聚類指標變化值?？梢钥闯觯財?shù)量為9 時，DBI為最低值，此時最小簇間距離第1 次明顯增大，而最大簇內(nèi)平均距離則變化不大，這意味著最后1 對間距較低的簇已完成合并；而當(dāng)簇數(shù)量減小至8 時，最大簇內(nèi)平均距離大幅升高，而最小簇間距離變化不大，說明分類數(shù)已經(jīng)過低，錯誤地將2 個距離較遠的簇合并在一起。綜上所述，DBI 最小時，意味著聚類進行到了1 個“簇內(nèi)相似度最高，簇間區(qū)分度最高”的最佳狀態(tài)，可以較好地表征任意測量參數(shù)組合下可區(qū)分的故障的數(shù)量。

圖1 3 項聚類度量指標與簇劃分數(shù)的關(guān)系Fig.1 Relationship between three clustering metrics and clustering number

2.2 魯棒性指標的計算

實際產(chǎn)品由于制造公差、環(huán)境因素、測試干擾等原因，遙測結(jié)果與仿真存在偏差，而DBI 只能衡量緊湊度和分離度的相對大小。如果傳感器排布方案對噪聲敏感度過高，會使最優(yōu)劃分數(shù)在實際工作環(huán)境中相對仿真模型得出的理論值發(fā)生偏移，導(dǎo)致評價依據(jù)失效。

噪聲對凝聚層次聚類的作用機理比較復(fù)雜，因此采用蒙特卡洛方法評估任意測量特征組合的抗噪性能。定義噪聲偏移率e為N次隨機噪聲試驗中，最優(yōu)劃分數(shù)發(fā)生變化的次數(shù)占總試驗次數(shù)的比率，其表達式為

式中：N為試驗次數(shù)；Ki為第i次試驗中，加入噪聲后的最優(yōu)劃分數(shù)；下標ζ代表加入隨機噪聲的均方根值。圖2 為2 組測點方案的魯棒性對比，其中RMS 代表噪聲的均方根值。圖2（a）表示2 組測點組合在不同強度隨機噪聲下進行100 次試驗后的最優(yōu)劃分數(shù)的偏移率。可以看出，組合2 的偏移率明顯低于組合1；從圖2（b）、圖2（c）可以看出，加入噪聲后，2 個組合的DBI 變化趨勢均有變化，但組合2 的DBI 曲線形狀基本保持不變，而組合1 的DBI 曲線在最小值附近發(fā)生了劇烈的變化。該結(jié)果說明，即使2 組排布對應(yīng)的最優(yōu)劃分數(shù)相近，但是其對噪聲的靈敏度可能有非常顯著的差別，這也表明了采用魯棒性指標的正確性和必要性。

圖2 2 組測點方案的魯棒性對比Fig.2 Robustness comparison of two measuring parameters schemes

綜上所述，結(jié)合多臺發(fā)動機的性能偏差先驗信息，本文選取100 次0.01 均方根噪聲下的噪聲偏移率e0.01作為評價指標。

2.3 可靠性指標的計算

飛行器對可靠性和安全性有著較高的需求，傳感器的布置應(yīng)該盡可能考慮其風(fēng)險系數(shù)。因此，需對發(fā)動機測量系統(tǒng)開展失效模式影響分析（FMEA），通過基于知識的方法確定傳感器組合的風(fēng)險程度。本文采取風(fēng)險指數(shù)矩陣作為風(fēng)險的量化指標，其可以作為故障嚴酷程度和發(fā)生概率的綜合評價指標。風(fēng)險指數(shù)矩陣的定義如表2所示。

表2 航天產(chǎn)品風(fēng)險指數(shù)矩陣Table 2 Risk index matrix of aerospace product

表3 發(fā)動機傳感器FMEA 表Table 3 FMEA of liquid rocket engine sensors

表4 不同權(quán)重參數(shù)的尋優(yōu)性能對比Table 4 Comparison of optimization performance between different weights

以發(fā)動機系統(tǒng)為初始約定層次，對測量參數(shù)對應(yīng)的傳感器的進行失效模式影響分析，確定故障模式、3 層次影響及風(fēng)險量化指標，其形式如表（3）所示。統(tǒng)計n種傳感器風(fēng)險指數(shù)，并形成n×1 的傳感器風(fēng)險指數(shù)向量Rn×1。

對于1 項發(fā)動機性能參數(shù)，其可能由單個傳感器直接測量，也可能需要多個傳感器經(jīng)過特定的解算才能獲得。定義1 個由0、1 組成的m×n的傳感器測量矩陣θm×n，其第i行、第j列元素代表第i種發(fā)動機性能參數(shù)是否需要安裝第j類傳感器?？紤]到泵壓式發(fā)動機內(nèi)組件無冗余備份，均為單點故障，因此可用單機級風(fēng)險指數(shù)之和作為該方案的總風(fēng)險量化指標，其表達式為

式中：Rtot為總風(fēng)險量化指標；θ為傳感器測量矩陣；R為傳感器風(fēng)險指數(shù)向量。

2.4 基于多目標二進制粒子群的測點組合優(yōu)化設(shè)計

由上述分析可知，測點組合的優(yōu)化目標為［-K（s），e（s），R（s）］。同時，為防止產(chǎn)生過多的Pareto 解，提高分析最優(yōu)解的效率，對3 項目標函數(shù)增加了約束。綜上，該優(yōu)化問題可表示為

式中：Kmin、emax和Rmax分別優(yōu)化過程中需滿足的最小劃分數(shù)、最大偏移率和最大風(fēng)險系數(shù)。為采用外點罰函數(shù)法將含約束問題轉(zhuǎn)換為無約束問題：

其中：v（s）為懲罰函數(shù)，其形式為

其中：G1～G3為每項測量性能指標的懲罰函數(shù)；κ1～κ3為懲罰函數(shù)。該問題屬于典型的組合優(yōu)化問題，求解難度較大，采用一般的方法難以得到準確的結(jié)果。本文采用改進的二進制多目標粒子群算法［22］優(yōu)化特征子集，根據(jù)實際需求人工篩選Pareto 前沿中的非支配解，確定最終的測量參數(shù)組合。其步驟如下：

步驟1初始化粒子群中的各粒子參數(shù)。隨機生成各粒子的位置S=｛s1，s2，…，sm｝和速度V=｛v1，v2，…，vm｝。

步驟2對各粒子按第3 步所述方法進行適應(yīng)度計算，按適應(yīng)度向量的支配關(guān)系確定個體最優(yōu)位置Pbest及全局最優(yōu)位置Gbest。

步驟3按式（30）和式（31）更新每個粒子的速度和位置。

式中：vi，k、xi，k為第i個粒子在第k次迭代中的速度向量和位置向量；c1、c2為認知系數(shù)；γ1、γ2、r為［0，1］之間的隨機數(shù)；ω為慣性權(quán)重。為進一步降低粒子群進入局部最優(yōu)解的風(fēng)險，引入Sigmoid 函數(shù)［23-24］對權(quán)重進行自適應(yīng)調(diào)整，其數(shù)學(xué)表達式為

其中：A、B為調(diào)節(jié)參數(shù)；ωmax、ωmin分別為慣性權(quán)重的最大、最小值；l為迭代次數(shù)；A、B為常數(shù)。

步驟4針對每個粒子，若更新后的粒子支配目前的個體最優(yōu)解Pbest，則更新Pbest；若與目前個體最優(yōu)解互不支配，則按50%概率更新Pbest。對Gbest的更新策略與Pbest相同。

步驟5將新的Gbest加入外部精英檔案并維護，去除檔案內(nèi)重復(fù)的粒子，并對支配關(guān)系進行排序，刪除全部被支配的粒子，確保檔案內(nèi)的粒子全部處于Pareto 前沿。

步驟6若迭代次數(shù)未達到200 次，重復(fù)步驟3～步驟5，迭代次數(shù)加1；否則，停止迭代，從外部檔案中取出全部解，轉(zhuǎn)化成測點排布方案，并根據(jù)設(shè)計需求，人工篩選測量方案。

綜上所述，基于多目標二進制粒子群算法的測點組合優(yōu)化方法流程圖如圖3 所示。

圖3 測量特征優(yōu)化流程圖Fig.3 Flow chart of measuring parameters optimization

3 算例分析

某現(xiàn)役、開式循環(huán)火箭發(fā)動機的測量參數(shù)為：燃料泵后壓力、氧化劑泵后壓力、渦輪轉(zhuǎn)速和氧化劑噴前壓力。該發(fā)動機的拓撲結(jié)構(gòu)如圖4所示。

圖4 某現(xiàn)役火箭發(fā)動機系統(tǒng)圖Fig.4 System diagram of rocket engine on active service

根據(jù)第1 節(jié)所述方法對發(fā)動機靜態(tài)特性建模并計算故障特征矩陣，獲得了L=510 種故障下m=37 項發(fā)動機性能參數(shù)的變化情況。采用2.1～2.3 節(jié)所述方法計算現(xiàn)有測點組合下性能指標，結(jié)果表明該組測點最優(yōu)劃分數(shù)為9，即可以分辨出9 種故障，噪聲偏移率為0。為更好地表征故障特征在空間中的分布，采用主成分分析（PCA）將故障特征降至三維。降維后的故障庫定義為

式中：變換矩陣Λ為歸一化的故障特征協(xié)方差矩陣Φ最大的3 個特征值對應(yīng)的特征向量構(gòu)成的矩陣。Φ的定義為

其中：μs為故障特征矩陣列向量的均值排列形成的行向量。降維后的故障空間特征如圖5 所示?？梢钥闯?，不同故障中心間隔距離較大，說明該排布的抗噪能力較強；但5 號故障中心聚集了271 種故障，且均與副系統(tǒng)相關(guān)。其原因是2 個泵后壓力測量值與渦輪轉(zhuǎn)速相關(guān)性很高，相當(dāng)于引入了冗余測點，測點的魯棒性有所提高，但無法區(qū)分副系統(tǒng)故障。

圖5 原測點布置方案下的故障特征主成分圖Fig.5 Principal components of fault characteristics under original measuring parameters scheme

圖6 為故障庫有效性的驗證過程。某次飛行故障中采集到的4 項測量參數(shù)的時域曲線如圖6（a）所示。取故障前、后的穩(wěn)態(tài)值，計算變化率并做規(guī)范化處理，計算故障特征到9 個故障分類簇中心的Euclid 距離，如圖6（b）所示；找到距離最近的1 個，即為最可能發(fā)生的故障類型。從結(jié)果可以看出，該故障被識別為副系統(tǒng)性能下降故障；通過遙測錄像等多媒體記錄進一步排查，確認為渦輪泵結(jié)構(gòu)破壞，為副系統(tǒng)性能下降故障的1 個子集。以上分析驗證了故障特征矩陣的準確性。

圖6 發(fā)動機靜態(tài)故障特征庫的驗證Fig.6 Verification of engine static fault feature list

采用第2 節(jié)中所述方法對傳感器排布進行優(yōu)化，其中，Kmin=8，Emax=30%，Rmax=0.6，粒子數(shù)量30，最大迭代數(shù)量100，ωmax=0.9，ωmin=0.4，A=13，B=0.15。圖7 為算法收斂過程，可以看到在44 次迭代之后，非支配解數(shù)量趨于穩(wěn)定。為驗證自適應(yīng)權(quán)重方法的有效性，選取30 組隨機初始粒子位置，對第1 個最優(yōu)非支配解出現(xiàn)的平均迭代次數(shù)I1和獲得全部最優(yōu)非支配解的平均迭代次數(shù)I2進行計算。其中，I1表征了算法尋優(yōu)的效率，I2則反映出算法的全局尋優(yōu)性能。計算結(jié)果如表（4）所示。結(jié)果表明，自適應(yīng)權(quán)重的I1整體上與ω=0.9 和ω=0.4 的非自適應(yīng)權(quán)重方法持平，有小幅度的降低；自適應(yīng)權(quán)重的I2則明顯優(yōu)于ω=0.9 的情形，略優(yōu)于ω=0.4 的情形。上述結(jié)果驗證了自適應(yīng)權(quán)重方法兼具搜索效率和全局尋優(yōu)性能。

圖7 多目標粒子群算法迭代過程Fig.7 Iteration process of MOBPSO

優(yōu)化算法收斂后共獲得8 組非支配解。8 組解非支配解中，最優(yōu)劃分數(shù)在9～13，偏移率均在20%以內(nèi)。8 組非支配解的空間分布情況和非支配解各性能參數(shù)分布情況如圖8、圖9 所示。

圖8 非支配解空間分布情況Fig.8 Spatial distribution of non-dominated solutions

圖9 非支配解各性能參數(shù)分布情況Fig.9 Distribution of performance parameters of nondominated solution

選取最優(yōu)劃分數(shù)13 的2 組解中的1 組進行對比分析，對應(yīng)的測點組合為：副系統(tǒng)混合比、主系統(tǒng)混合比、轉(zhuǎn)速和氧噴前壓力。該選取采用主成分分析對故障特征進行降維，3 項主成分的分布圖如圖10 所示?？梢钥闯觯瑢y點排布進行優(yōu)化后，可分辨的故障從9 種增加至13 種；各故障中心兩兩之間距離較遠，偏移率約為3.1%，與原排布方案相當(dāng)。其中，副系統(tǒng)故障的分辨能力明顯提升，副系統(tǒng)故障中心（7 號）雖仍存在一定程度的聚集，但氧、燃副路堵塞和氧副路泄漏這3 種重要故障已經(jīng)從中剝離出來。上述結(jié)果表明了優(yōu)化方法的有效性。

圖10 優(yōu)化布置方案下的故障特征主成分圖Fig.10 Principal components of fault characteristics under measuring parameters scheme

進一步觀察優(yōu)化過程，發(fā)現(xiàn)包含副系統(tǒng)混合比的測量方案的可分辨故障數(shù)普遍高于不包含副系統(tǒng)混合比的測量方案。因此采用蒙特卡洛方法，隨機選擇200 組包含和不包含副系統(tǒng)混合比參數(shù)的4 測點方案并統(tǒng)計其3 項指標，結(jié)果如圖11 所示?？梢钥闯觯毕到y(tǒng)混合比的測量方案的最優(yōu)劃分數(shù)偏高，偏移率變化不大，風(fēng)險指數(shù)略有上升。圖中紅色部分代表包含“副系統(tǒng)混合比”的測量方案；藍色部分代表代表不包含“副系統(tǒng)混合比”的測量方案。

圖11 副系統(tǒng)混合比對故障分辨能力的影響Fig.11 Influence of mixing ratio of subsystems on fault recognition

從發(fā)動機機理上可以解釋上述現(xiàn)象產(chǎn)生的原因。作為副系統(tǒng)氧、燃流量的比值，副系統(tǒng)混合比直接表征了副系統(tǒng)的工作狀態(tài)；另外，從機理上來看，副系統(tǒng)混合比直接決定了燃氣發(fā)生器內(nèi)高溫氣體的組分、溫度和壓力，進而顯著影響渦輪入口燃氣焓流量和渦輪泵輸入功率；而在發(fā)動機主系統(tǒng)中，渦輪泵的揚程一旦發(fā)生改變，壓力平衡方程的解也會發(fā)生變化，從而導(dǎo)致發(fā)動機靜態(tài)參數(shù)相較額定工況發(fā)生偏移。由于副系統(tǒng)管徑小、壓力高，目前副系統(tǒng)流量一般采用軟測量方法。由此可見，持續(xù)開展高精度的副系統(tǒng)流量測量硬件的研制和測量方法的研究是改善靜態(tài)故障的辨識率的重要方向。

綜上所述，所提方法能夠提高測量系統(tǒng)對故障的分辨能力，且兼具良好的魯棒性和系統(tǒng)可靠性。后續(xù)工程化過程中，主要還需要考慮硬件、軟件兩方面。對于硬件，尤其是參與閉環(huán)控制的傳感器，一旦發(fā)生故障可能導(dǎo)致發(fā)動機虛警、執(zhí)行錯誤動作甚至觸發(fā)備保關(guān)機，屬于Ⅱ類甚至Ⅰ類單點故障。因此，應(yīng)采取冗余措施，例如多傳感器冗余表決或根據(jù)FMEA 在控制器模塊制定IF-THEN 形式的傳感器失效后處置策略；此外，對傳感器測量值采取合適的濾波方法，也是在線故障識別的重點、難點之一。對于軟件和算法層面，主要需要考慮軟測量下的誤差量化。對于壓力、轉(zhuǎn)速等信息，可以直接采用傳感器測量；而本文所提提到的“副系統(tǒng)混合比”指標，則需要副系統(tǒng)管路某節(jié)流元件兩端壓力和水試阻力系數(shù)實時解算流量，再計算混合比。水試阻力系數(shù)測試不合理或未覆蓋飛行中的壓差將會顯著影響故障診斷的可靠性。因此，應(yīng)開展冷態(tài)流動試驗，對副系統(tǒng)管路和節(jié)流元件的流阻特性進行精確的標定，且工況覆蓋飛行中可能出現(xiàn)的溫度、壓力。目前，火箭發(fā)動機在線故障判別方法已在現(xiàn)役常溫運載火箭實現(xiàn)了初步應(yīng)用。

4 結(jié)論

1）基于火箭發(fā)動機靜態(tài)特性的數(shù)學(xué)模型，通過輸入模擬故障獲得發(fā)動機靜態(tài)特性故障庫。依托基于模型的故障特征庫，考慮故障分辨率、測量的魯棒性和系統(tǒng)可靠性，提出了液體火箭發(fā)動機測量特征子集評價指標的計算方法以及基于改進的多目標二進制粒子群算法的傳感器測點組合優(yōu)化方法。

2）將上述方法應(yīng)用于某現(xiàn)役火箭發(fā)動機的測量特征優(yōu)化。優(yōu)化后的測點布置方案，其可分辨的故障數(shù)量從9 種提升至13 種，魯棒性變化較小，風(fēng)險指數(shù)有較小提升；進一步探究了“副系統(tǒng)混合比”這一測量參數(shù)在泵壓式發(fā)動機故障分辨中的重要作用，并從機理層面進行解釋。本文提出的方法對其他復(fù)雜、閉環(huán)動力系統(tǒng)測量特征的選擇具有較好的應(yīng)用價值。