基于動態故障樹的衛星可靠性分析

王許煜，胡敏，張學陽，趙玉龍，李玖陽

航天工程大學，北京 101400

隨著國內外航天技術的快速發展，衛星已經被廣泛應用于通信遙感、導航定位和軍事偵察等領域。然而作為一類特殊產品，衛星具有其他產品不具備的特點[1-2]：一是衛星設計和制造周期長，成本高且生產數量少；二是衛星需要苛刻的發射條件，且會長期受到復雜惡劣的太空環境的影響；三是若衛星在軌運行期間發生故障，無法對其進行有效的人為維修。因此，通常把衛星稱為四高產品，即高技術、高投入、高效益和高風險產品[3]。對于衛星星座而言，衛星可靠性關系到星座可用性和備份策略的選取[4]，所以為了能夠連續穩定地完成特定任務，需要對衛星可靠性進行分析，以確保衛星的高可靠性。

目前傳統的衛星可靠性模型通常采用指數模型，然而該簡化模型雖然有利于分析問題，但不能準確地描述衛星早期故障和損耗累積的過程。文獻[5-6]利用多個概率分布的組合來表示衛星各個壽命階段進行可靠性建模。文獻[7]借鑒GPS衛星的經驗，分別采用Weibull分布模型和正態分布模型模擬衛星隨機故障和損耗故障，進而建立衛星可靠性模型。利用分布模型建模雖然能反映衛星的實際情況，但由于衛星系統可靠性高，失效樣本十分有限，導致模型參數難以確定，不利于可靠性模型的準確建立。文獻[8]提出基于序貫優化和可靠性分析的交叉熵優化方法，對衛星可靠性設計進行優化分析。文獻[9]利用貝葉斯網絡模型提出一種基于不確定隸屬度函數和區間特征量的復雜不確定系統可靠性分析方法。這些方法雖然能有效分析衛星可靠性，但都無法準確描述衛星系統內部各組件之間復雜的動態邏輯關系。本文利用動態故障樹方法對衛星系統進行分析，為其可靠性設計分析提供了一個新的思路。

傳統的故障樹分析法(fault tree analysis，FTA)是一種分析系統可靠性的有效方法[10]，但其不適用于具有動態失效行為的復雜系統。馬爾可夫模型雖然可以表達系統的動態特性，但是由于模型的狀態空間數量會隨系統規模的增長呈指數增長，而導致模型的建立和求解十分復雜[3]。動態故障樹(dynamic fault tree，DFT)分析方法結合了傳統FTA分析方法和Markov模型兩者的優點，它通過引入描述系統動態特征的邏輯門建立相應的動態故障樹，并將其轉化為馬爾可夫模型，從而能夠對動態系統進行可靠性和安全性分析，被廣泛應用于動態系統的可靠性建模[11]。

本文首先采用馬爾可夫鏈和二元決策圖相結合的分析方法，建立衛星的電源、姿軌控和推進3個分系統的動態故障樹模型，在此基礎上得到衛星的隨機故障模型，并綜合考慮損耗故障建立衛星可靠性模型。接著利用蒙特卡洛仿真對隨機故障模型進行評估分析，最后將其與Weibull分布模型進行性能比較。

1 動態故障樹理論

1.1 動態故障樹

故障樹是一種特殊的樹狀邏輯因果關系圖，它用規定的事件、邏輯門和其他符號描述系統中各事件之間的因果關系[12]。動態故障樹是指包含一個或多個動態邏輯門的故障樹，與傳統的故障樹相比，由于其引入了能夠表征系統動態時序特性的元素，可以用于復雜冗余系統的分析求解。

動態故障樹中的動態邏輯門是描述事件之間因果關系的邏輯符號，其能夠用來表征與事件發生順序相關的動態特性，這些特性必須考慮不同底事件的發生順序以及各組件之間的依賴關系[12]。動態邏輯門主要有：優先與門、順序相關門、功能相關門、冷備件門、溫備件門和熱備件門等[13]。

(1)優先與門

優先與門(priority-AND gate, PAND)具有兩個底事件A和B，只有這兩個底事件按照從左到右的順序依次發生時才會導致頂事件的發生，圖1表示優先與門轉化為馬爾可夫鏈的過程。其中，狀態的第一個數字表示底事件A的狀態，第二個數字表示底事件B的狀態，數字1和0分別表示底事件故障和正常的兩種狀態。Fa和Op分別表示頂事件T發生和不發生。轉移路徑上的符號表示相應事件發生故障，其轉移率為該事件的失效率。

(2)順序相關門

順序相關門(sequence enforcing gate, SEQ)具有多個底事件，當底事件按照從左到右的順序依次發生時頂事件T才發生。以3個底事件為例說明順序相關門轉化為馬爾可夫鏈的過程，如圖2所示。

圖2 順序相關門轉化為馬爾可夫鏈的過程Fig.2 The process of transforming sequence enforcing gate into Markov chain

(3)功能相關門

功能觸發門(functional dependency gate, FDEP)表示系統中各個事件之間的關系，并沒有實際的輸出，若系統中的觸發事件T發生，則與其相關的所有基本事件A和B都發生，而基本事件的發生對觸發事件的發生則沒有任何影響。圖3為功能相關門轉化為馬爾可夫鏈的過程。

圖3 功能相關門轉化為馬爾可夫鏈的過程Fig.3 The process of transforming functional dependency gate into Markov chain

(4)冷備份門

冷備件門(cold spare gate, CSP)包括一個主件和至少一個冷備件。系統運行時主件處于正常工作狀態，其余的冷備件處于不工作狀態，且失效率為零。當主件A發生故障時，冷備件開始工作，直至所有冷備件失效，頂事件T才會發生。圖4冷備份門轉化為馬爾可夫鏈的過程。

圖4 冷備件門轉化為馬爾可夫鏈的過程Fig.4 The process of transforming cold spare gate into Markov chain

(5)溫備份門

溫備件門(warm spare gate, WSP)的邏輯關系與冷備件門基本相同，不同的是當系統主件工作時，溫備件處于預工作狀態具有一定的失效率，但小于正常工作狀態下的失效率, 只有主件與溫備件都發生故障時頂事件T才會發生。圖5是溫備件門轉化為馬爾可夫鏈的過程。

圖5 溫備件門轉化為馬爾可夫鏈的過程Fig.5 The process of transforming warm spare gate into Markov chain

(6)熱備份門

熱備件門(hot spare gate, HSP)中的主件A和熱備件都處于正常工作狀態，熱備件具有正常工作狀態下的失效率。當主件與熱備件都發生故障時頂事件T才會發生。圖6是熱備件門轉化為馬爾可夫鏈的過程。

圖6 熱備件門轉化為馬爾可夫鏈的過程Fig.6 The process of transforming hot spare gate into Markov chain

1.2 動態故障樹的分析方法

在動態故障樹模型中通常只有很少的邏輯門具有動態特性。所以在動態故障樹分析中，首先利用線性搜索算法對DFT進行模塊化分解，該算法通過對DFT進行兩次深度優先最左搜索(DFLM)[12]，得到相互獨立的動態子樹和靜態子樹，其中包含動態邏輯門的子樹為動態子樹，只有靜態邏輯門的子樹則為靜態子樹。通過模塊化處理后，復雜系統的DFT被分解成多個子模型，極大簡化了模型的求解過程[14]。然后對不同類型的子樹單獨求解，利用二元決策圖(binary decision diagram，BDD) 對靜態子樹進行求解, 而動態子樹則轉化成相應的Markov 鏈進行求解。最后結合各獨立子樹分析求得的故障概率進行綜合分析，通過邏輯關系自下而上分步求解直至獲得動態故障樹的頂事件發生概率。動態故障樹的整個分析過程如圖7所示。

圖7 動態故障樹模型分析過程Fig.7 Dynamic fault tree model analysis process

2 衛星可靠性模型分析

衛星主要由航天器平臺和有效載荷兩大部分組成，其中航天器平臺包括結構分系統、電源分系統、數管分系統、測控分系統、姿軌控分系統、推進分系統以及熱控分系統。圖8為衛星的簡要故障樹。根據衛星在軌故障的統計結果顯示[15-16]，對于衛星而言，電源、姿軌控和推進3個分系統在軌故障所占比例最高，而且一旦該系統發生故障將直接影響航天器的正常運行。本文重點針對這3個分系統展開分析。

圖8 衛星簡要故障樹Fig.8 Satellite fault tree

2.1 動態故障樹建模

(1)電源分系統

目前，衛星上大多采用太陽電池翼-蓄電池組的電源系統，由兩條相同的供電母線并聯組成，以提高系統的可靠性，每條母線包括蓄電池組、太陽電池陣以及電源控制器。其中，蓄電池組采用鋰離子蓄電池組，該蓄電池組由4組蓄電池并聯組成，每組蓄電池由9個單體串聯而成，并假設每個單體都一樣。電源系統的動態故障樹模型如圖9所示。

圖9 電源系統的動態故障樹模型Fig.9 Dynamic fault tree model of power system

模型中各事件說明如表1所示。

表1 電源系統動態故障樹模型事件說明

(2)姿軌控分系統

衛星姿軌控分系統多數采用精度高、穩定性好的三軸姿態穩定控制方式，主要由控制器子系統、執行機構子系統和姿態確定子系統組成，它們之間通過相互作用實現衛星姿態的閉環控制。其中，執行機構采用三正一斜裝的反作用輪，并利用磁力矩器和推力器對其進行卸載。姿態確定系統則采用陀螺和星敏感器的組合，同時將太陽敏感器和兩個紅外地平儀作為星敏感器的備份。姿軌控系統的動態故障樹模型如所圖10所示。

模型中各事件說明如表2所示。

(3)推進分系統

衛星采用雙組元統一推進系統，其技術成熟，被廣泛應用于長壽命航天器上。該系統利用氣瓶將燃燒劑(MMH)和氧化劑(MON-1)以霧化的形式推出，并通過氣管和各類閥門進入燃燒室，在按一定的比例燃燒后產生推力。推進分系統的故障樹模型如圖11所示。

模型中各事件說明如表3所示。

圖10 姿軌控分系統的動態故障樹模型Fig.10 Dynamic fault tree model of attitude and orbit control system

表2 推進分系統動態故障樹模型事件說明

圖11 推進分系統的故障樹模型Fig.11 Fault tree model of propulsion system

表3 推進分系統故障樹模型事件說明

2.2 基于BDD的靜態子樹分析

衛星系統動態故障樹模型規模較大，本文將以電源系統為例詳細分析動態故障樹模型的求解過程，其余兩個系統將以同樣的方法進行分析。由于兩條供電母線相互獨立且組成相同，所以可對供電母線1進行單獨分析，將動態故障樹進行模塊化分解，得到相應的靜態子樹和動態子樹，其中靜態子樹包括G31、G41和G51，動態子樹包括G131。

對于靜態子樹，采用BDD進行求解。BDD 是一個有向無環圖[17]，由根節點、葉節點和中間節點組成。每條BDD路徑從根節點開始，經過中間節點時進行布爾函數的計算，并有兩條分支，分別代表事件失效和正常兩種狀態，最終到達葉節點。葉節點有兩種類型，其中“1”表示系統故障，“0”表示系統正常。

(1)G31靜態子樹

將G31靜態子樹進行BDD分解首先要確定底事件的最優指標順序，為便于計算，假設指標順序為：E11

根據BDD可以得G31事件發生的概率為:

P(G31)=P(E11∪E21∪E31∪E41∪E51∪

E61∪E71∪E81∪E91)=1-

(1)

圖12 G31靜態子樹的BDDFig.12 BDD for G31 static subtree

(2)G41靜態子樹

由于G41靜態子樹由4個相同的中間事件并聯所得，所以為了簡化分析，先對其中一個中間事件G111進行單獨分析，然后再結合其余事件進行綜合分析。G111由9個底事件串聯所得，所以可以根據其組成結構得到G111發生的故障概率為：

P(G111)=1-[1-P(E101)]9

(2)

然后將G91、G101、G111、G121當作G41的底事件進行分析，假設底事件指標順序為：G91

圖13 G41靜態子樹的BDDFig.13 BDD for G41 static subtree

根據BDD可以得G41事件發生的概率為:

P(G41)=P(G91G101G111G121)

(3)

(3)G51靜態子樹

對于G51靜態子樹，由于熱備份件門的存在，首先要將動態子樹G131作為其輸入底事件，假設底事件指標順序為：E191

根據BDD可以得G51事件發生的概率為:

P(G51)=P(E191∪E201∪G131)=

(4)

最終，結合各模塊事件的發生概率，可以得到供電母線1的故障概率，并進一步得到電源系統T11的故障概率。

P(T11)=P(G11G21)=P(G11)P(G21)=

(5)

2.3 基于Markov的動態子樹分析

根據動態故障樹模塊化分析可得，電源動態故障樹中只有G131為動態子樹，對于動態子樹，將采用Markov鏈進行求解。Markov鏈是用來表示系統各狀態之間轉換的一種方法，由于系統下一時刻的狀態只與當前時刻的狀態有關，而與之前任何時刻無關，所以具有“無記憶性”。在Markov 鏈中圓圈代表系統當前時刻的狀態, 有向線段表示狀態的轉移方向，而上面的權值則表示狀態的轉移概率。

求解馬爾可夫鏈的方法有很多，文獻[18]提出將系統的Markov狀態轉移圖分解成多條狀態轉移鏈，并推導出不同鏈長的轉移鏈的概率計算公式，進而求得動態子樹的發生概率。鏈長為n的狀態轉移鏈Tn的概率計算公式：

(6)

式中：λi-1,i>0，為狀態i-1到狀態i的轉移率；λi-1,Op≥0，為狀態i-1到狀態Op的轉移率。

由圖6可得熱備件門具有兩條鏈長為2的狀態轉移鏈，根據式(6)計算各鏈的概率，然后將其相加，就可得到動態子樹G131的發生概率，最后將所得概率代入(4)中進而得到電源系統的故障概率。

2.4 結果分析

假設電源系統各部件的壽命均服從指數分布，且各底事件的失效率如表4所示。設衛星在軌工作時間為105h，壽命末期為8.76×104h。

表4 電源分系統各底事件失效率

將相關數據代入模型，得到衛星電源分系統的可靠性隨時間變化曲線，如圖15所示。

圖15 電源分系統的可靠性隨時間變化曲線Fig.15 Reliability curve of power system with time

利用動態故障樹方法分別對姿軌控分系統和推進分系統進行分析，并最終得到衛星隨機故障模型的可靠性隨時間變化曲線，如圖 16所示。由圖 16可知，隨著衛星工作時間的增長，衛星及3個分系統的可靠性都在逐漸下降，同時電源分系統對衛星可靠性的影響在逐漸增大，到衛星壽命末期時，成為3個分系統中最為薄弱的環節。

圖16 衛星隨機故障模型的可靠性隨時間變化曲線Fig.16 Reliability curve of satellite random fault model with time

衛星隨機故障模型在不同時間的可靠性如表5所示。利用動態故障樹對衛星分系統進行建模分析，能夠有效地得到各系統以及整星的隨機故障發生概率，從而便于分析系統的薄弱環節，進一步提高衛星的可靠性。

最后，綜合考慮衛星的損耗故障，將隨機故障模型與損耗故障模型相乘得到衛星的可靠性模型，本文利用正態分布模型模擬衛星損耗故障[19]，根據衛星壽命末期的要求，確定模型的期望為93 000，標準差為10 000，得到衛星的可靠性模型如圖17所示。由圖17可知，在衛星早期階段，隨機故障對衛星可靠性的影響較大，衛星可靠性隨著時間的增大在緩慢下降，當衛星工作6×104h后，可靠性開始迅速下降，此時損耗累積的作用逐步顯著，到達衛星壽命末期8.76×104h時，衛星的可靠性為0.59，此時損耗累積已成為影響衛星可靠性的主要因素，至工作時間105h時，衛星可靠性僅為0.2。

表5 衛星隨機故障模型在不同時間的可靠性

圖17 衛星可靠性模型Fig.17 Satellite reliability model

3 基于蒙特卡洛仿真的模型評估

蒙特卡洛仿真是分析動態故障樹的重要方法，其具有較高的計算精度，但計算時間較長，同時該方法也是評估其他方法的標桿。本文利用蒙特卡洛仿真方法對所得的衛星隨機故障模型進行評估分析。仿真時，根據表4中底事件服從的壽命分布及參數，對各事件出現的時間進行抽樣，將所得樣本作為仿真時鐘的推進點，并在推進點對系統進行評估和參數計算，直至系統失效或到達最大仿真時間，在仿真次數為105的條件下得到相應的可靠性曲線如圖 18所示。

圖18 基于蒙特卡洛仿真的可靠性曲線Fig.18 Reliability curve based on Monte Carlo simulation

在不同時間下，將本文方法與基于蒙特卡洛仿真方法進行對比，得到的計算結果如表6所示。結果表明，本文采取的動態故障樹分析方法能夠有效的分析衛星的隨機故障，且所得的隨機故障模型具有良好的精確性。

表6 本文方法與基于蒙特卡羅仿真方法對比

4 Weibull分布模型比較

4.1 Weibull分布參數擬合

Weibull分布模型具有較強的適應性，能夠擬合不同類型的數據，被廣泛應用于系統可靠性的分析。本文在基于蒙特卡洛仿真所得到的可靠性數據的基礎上，利用最大化擬合優度參數估計方法擬合得到Weibull分布模型。擬合優度表示擬合數據與模型之間的差異，其值越大則表明模型擬合度越好。基于蒙特卡洛仿真得到的可靠性擬合數據記為R(ti)，Weibull分布模型的可靠性記為R′(ti,α,β)，其中α為尺度參數，β為形狀參數，則擬合優度R2可通過以下公式得到：

(7)

(8)

(9)

(10)

利用遺傳算法對擬合優度函數進行優化，得到擬合優度的最大值，進而確定Weibull分布模型的參數。最終Weibull分布模型尺度參數α為295 742，形狀參數β為1.46，擬合優度R2為0.999 3。

4.2 模型比較

以基于蒙特卡洛仿真所得的可靠性數據為依據，將擬合得到的Weibull分布模型與本文所得的隨機故障模型進行對比，計算兩種模型與可靠性數據之間的殘差絕對值,如圖 19所示。由圖 19可知，Weibull分布模型的殘差絕對值的峰值基本上都大于本文所得的隨機故障模型，并進一步求得殘差絕對值的平均值，分別為0.001 2和0.000 7。

圖19 兩種模型與可靠性數據之間的殘差絕對值Fig.19 Absolute value of residual between two models and reliability data

結果表明，基于動態故障樹方法所得的隨機故障模型能更好地貼合可靠性數據，更有效地模擬衛星部件隨機故障的動態行為。

5 結束語

本文采用馬爾可夫鏈和二元決策圖相結合的動態故障樹分析方法對衛星及3個分系統的動態邏輯和故障機制進行研究，建立完整的衛星動態故障樹，得到衛星隨機故障模型，并將隨機故障模型和正態分布模型結合得到更為準確的衛星可靠性模型，解決了傳統故障樹分析方法無法描述衛星系統動態特性的問題。同時，該方法有利于分析系統的薄弱環節，進而提高衛星的可靠性。而后在蒙特卡洛仿真數據基礎上對隨機故障模型和Weibull分布模型進行評估比較，結果表明本文方法能更準確地反映衛星系統隨機故障的動態行為，更利于衛星隨機故障模型的建立和分析。