蒙特卡洛法在發動機推力測量不確定度評估中的應用

鄭 科, 耿衛國,朱子環

(北京航天試驗技術研究所，北京 100074)

0 引言

發動機常被用于衛星、火箭等的精確軌道控制和姿態調整，姿軌控發動機推力矢量直接關系到衛星能否入軌以及發射任務的成??；準確測出推力矢量參數，能夠為發動機的在軌工作狀態提供基本依據。目前推力矢量測量有多種方案，但在我國的發動機高空模擬熱標定試驗中，技術和工藝比較成熟且經多次飛行驗證效果顯著的，是北京航天試驗技術研究所投入使用的轉臺推力矢量測量方法[1]和該所研制的基于壓電的動態矢量推力測量方法[2-3]。

無論采用何種推力測量方法，推力矢量參數都不是直接測得，而是基于一定的數學模型通過計算間接得出，前期對推力矢量不確定度的評估，是依據國家計量技術規范《測量不確定度評定與表示》JJF 1059.1-2012中的GUM方法，對各個具體測得量進行評估后，按照不確定度傳播率計算合成標準不確定度[4-5]。但針對多個輸入量和單一輸出量的測量模型，作為JJF 1059.1-2012的補充文件，國家計量技術規范《用蒙特卡洛法評定測量不確定度》JJF 1059.2-2012規范了一種評估測量不確定度的方法[6]，其核心內容是在建立測量模型的基礎上利用對概率分布的隨機抽樣進行分布傳播不確定度。本文基于壓電式推力矢量架，在介紹其數學模型后，采用GUM法和蒙特卡洛法(MCM,Monte Carlo method)進行不確定度評估，對不確定度評估結果進行對比驗證，比較兩種方法的適用性。

1 壓電式推力矢量測量

火箭發動機壓電式推力矢量系統的數學測量模型如圖1所示[7]。

圖1 某壓電式推力矢量測量系統數學模型

圖1中，O-XYZ為測力平臺坐標系，O′-XY′Z′為火箭發動機坐標系。其中，YOZ為壓電測力儀三維測力傳感器的安裝定位平面，O點為OX軸與三維測力傳感器安裝定位平面交點；Y′O′Z′為火箭發動機的安裝定位法蘭平面，O′點為火箭發動機噴管幾何理論軸線與法蘭平面交點。矢量力偏斜角為α，作用點為(δy,δz)。

F為姿軌控發動機點火產生的空間矢量推力；

Fx，Fy，Fz分別為矢量推力F在作用點(δy,δz)的三向分力；

Fx1，Fx2，Fx3，Fx4分別為通過測力平臺上的4個排列成正方形分布的三維力傳感器實際測得的X方向的4個分力。同樣，Fy1，Fy2，Fy3，Fy4和Fz1，Fz2，Fz3，Fz4分別為通過傳感器實際測得的Y方向和Z方向分力；

Mo為火箭發動機矢量推力F對測量平臺坐標系中心產生的總力距；

Mx，My，Mz分別為總力矩Mo在各方向上的分力矩；

a為發動機測力平臺上三維力傳感器與坐標軸之間的距離；

c發動機測力平面上三維力傳感器與法蘭平面之間的垂直距離。

建立火箭發動機推力矢量參數的測量模型：

發動機的三項推力為:

(1)

推力矢量力矩為:

(2)

則得:

推力斜偏角α:

(3)

側向力方位角:

(4)

推力偏移:

(5)

推力偏移方位角:

(6)

2 蒙特卡洛法不確定度評估

2.1 蒙特卡洛法原理

蒙特卡洛法計算測量不確定度是以概率統計為基礎，利用隨機抽樣實現概率分布傳遞的一種數值計算方法，區別于GUM法的不確定度傳播率，適用于可由任意多個概率密度函數(PDF,probability density function)表征的輸入量和單一輸出量的模型，尤其適用于各不確定度分量的大小不相近、輸出量的估計值和其標準不確定度相當、測量模型明顯呈非線性等典型情況[8-9]。

蒙特卡洛法為輸出量的表征提供了一種方法，輸出量Y的分布函數為。

(7)

式中，GY(η)為輸出量Y的概率密度函數。

gY(η)可由輸入量Xi的概率密度函數gxi(ξi)通過測量模型傳遞得到，概率分布傳遞如圖2所示，核心是對輸入量的PDF重復采樣，即利用對輸入量概率分布的隨機抽樣代入數學模型進行分布傳遞，計算求得輸出量Y的PDF的離散抽樣值，因為離散抽樣值GY(η)包括了輸出量Y的數值特性，所以能夠由輸出量的離散分布數值求出Y的最佳估計值以及標準不確定度，Y的不確定度和包含區間等數值特性的可靠性可隨著對輸入量的概率密度函數的隨機抽樣數M增加而提高[10-11]。

圖2 由多輸入量PDF得到單一輸出量PDF的示意

MCM評估過程：

1)建立輸出量Y和各個輸入量X1,X2,…,Xn的模型Y=f(X1,X2,…,Xn)；

2)確定利用獲得的先驗信息，確定各輸入量的概率密度函數gxi(ξi)；

3)確定MCM的仿真次數M；

4)從各個輸入量Xi的概率密度函數中隨機抽樣M個樣本值xir(i=1,2,…,N，r=1,2,…,M)；

5)對每個的樣本矢量x1r,x2r,…,xNr)，代入模型得yr=f(x1r,x2r,x3r,…,xNr)(r=1,2,…,M)；

6)將計算得到的M個模型值按照遞增順序排序，得出輸出量Y的離散表示G；

7)由輸出量的離散表示G計算Y的期望值以及標準差，求出估計值y和不確定度μ(y)；

8)通過輸出量的離散表示G求出給定包含概率下的包含區間[12]。

蒙特卡洛法的流程如圖3所示。

圖3 蒙特卡洛法流程圖

2.2 蒙特卡洛法的實現

2.2.1 概率密度函數的確定

針對輸入量數據量小的問題，通常采用最大熵原理[13-14]求出各輸入量的PDF，通過對每個輸入量的概率密度函數進行M次抽樣即解決了小樣本數據量小的問題。

2.2.2 抽樣仿真次數M

抽樣模擬仿真次數M增加，樣本容量的大小增加，蒙特卡洛法接近于輸出量的真實總體。但如果M增多，隨機抽樣需要計算的時間越久，M減小，與輸出量的實際情況偏離，計算不確定度的結果不準確。所以為了計算相對準確的不確定度需要合理的選擇仿真次數M。

抽樣仿真次數M的確定一般有兩種方法：

1)一般情況下，M的值應至少大于1/(1-p)的104倍，包含概率p為0.95時，抽樣仿真次數M:

(8)

2)一般要求測量結果的測量不確定度不超過兩位有效數字時:

(9)

根據自由度的計算公式可得:

(10)

因此，抽樣仿真次數M應該大于8×104，才能使評估結果的不確定度具有兩位有效數字。

根據上述兩種方法得到的計算結果，結合不確定度評估的實際經驗，一般可以取抽樣仿真次數M為106。

2.2.3 舍選法抽樣

在使用蒙特卡洛法評估測量不確定度時，需要在各輸入量PDF的約束下產生大量的隨機數，對測量模型進行隨機抽樣，即產生隨機變量，通過對產生的隨機數進行模型傳遞以及統計計算，求得輸出量的統計性質。

針對本文的小樣本數據，在通過最大熵原理求出針對小樣本數據的概率密度函數后，求得的概率密度函數分布不一定是常見分布，不能直接通過Matlab特定函數產生隨機數，為了得到符合任意概率分布下的隨機數，本文采用舍選抽樣法[15-16]求隨機數。

舍選抽樣法的原理是：根據給定的輸入量的概率密度函數f(x)，對概率分布為均勻分布的隨機數序列{R}進行舍選。舍選的原則是當f(x)取值較大時，選擇較多的隨機數ri；在f(x)取值較小時，選擇較少的隨機數ri，從而使最后得到的子樣本分布滿足給定的概率密度函數。

輸入量xi的概率密度函數為gxi(ξi)，在gxi(ξi)的約束下產生M個值xir(i=1,2,3,…,n,r=1,2,3,…,M)，若對n個輸入量分別在其概率密度函數的約束下產生M個偽隨機數，則可以得到M組向量:

(11)

2.2.4 蒙特卡洛法仿真結果

(12)

(13)

3 壓電式推力矢量不確定度評估

3.1 測量數據

根據某次推力試車數據推導得出推力矢量參數如表1所示。

表1 某次試車推力矢量參數

表2～4為試車前X、Y、Z方向的推力檢驗值。

表2 壓電式推力矢量系統X向靜態標定

3.2 GUM法不確定的評估

計算參數合成不確定度一般表達式為:

(14)

式中,f為被測量y與直接測得量xi的函數關系。

由于Fx=Fx1+Fx2+Fx3+Fx4，是線性的，所以對Fx的A類不確定度計算如下：

根據合成標準不確定度推導：

u2(Fx1) +u2(Fx2) +u2(Fx3) +u2(Fx4)

(15)

主推力合力的A類不確定度：

(16)

由前面章節推導公式得出推力矢量參數的合成不確定度。例如:

1)α的合成不確定度:

(17)

2)δ的合成不確定度:

(18)

通過上述表得出X、Y、Z方向的不確定度Ux，Uy，Uz，代入合成不確定度評估公式得出壓電式推力矢量參數(α,γ,β,δ)的不確定度評估結果如表5所示。

表5 GUM法推力矢量不確定度評估

表5中，擴展不確定度包含因子為k=2，包含區間的包含概率為95%。

3.2.1 GUM法的適用條件

在實際應用中，GUM法適用于以下條件：

1)輸入量的概率分布可以近似是為對稱分布；

2)輸出量的概率分布可以近似是正態分布或t分布；

3)測量數學模型為線性模型或者能夠用線性模型近似的模型[17]。

因此，GUM法具有一定的局限性。

當測量數學模型是復雜的非線性模型時，輸入量的一階偏導數通常不容易求解；根據泰勒級數展開，將非線性模型近似轉化為線性模型，忽略了高階級數項，帶來了一定的誤差，且每個輸入量之間的相關系數難以計算；評估擴展不確定度時，一般取包含因子k為2或3，具有主觀性，GUM法默認被測量的概率分布近似是正態分布或者t分布。

3.3 MCM的軟件實現

圖形用戶界面(GUI，graphical user interfaces)是Matlab的一個主要功能，是根據用戶體驗和用戶需求設計的可視化用戶界面，一般由窗口、菜單、對話框等各種圖形對象組成，用于與計算機、操作系統和應用程序互動交流，即使計算機產生某種動作，實現計算和繪圖等功能[18-19]。

本文根據上文中的蒙特卡洛法流程利用GUI編寫軟件，固化流程，通過輸入數據，進行直接計算，軟件不僅可以計算推力矢量參數的不確定度，通過更改數學模型Y=f(X1,X2,…,Xn)，還可以計算發動機試驗中溫度、流量以及壓力等的測量不確定度。

軟件具體實施方式：

1)將發動機試驗數據導入軟件界面；

2)通過多組數據計算每個輸入量的中心矩mi；

3)確定每個輸入量的數據邊界；

4)依據最大熵原理計算每組數據的概率密度函數gxi(ξi)，得出每個輸入量的概率分布；

5)利用舍選抽樣法進行抽樣，隨機抽樣符合各個輸入量Xi概率密度函數的M個樣本值(xir(i=1,2,…，6，r=1,2,…，M)；

6)點擊各輸出量即參數測量模型，得出輸出量的概率分布，從而給出不確定度及其給定包含概率下的包含區間。

將表3X方向的靜態標定線性差值，以及表4Y方向、表5Z方向的靜態標定線性差值導入不確定度評估軟件中，如圖4所示，不確定度評估結果如表6所示。

表3 壓電式推力矢量系統Y向靜態標定

表4 壓電式推力矢量系統Z向靜態標定

圖4 MCM推力矢量不確定度評估

表6 基于MCM的推力矢量不確定度評估結果

表6中，給出的包含區間為包含概率為95%的包含區間。

圖4中，導入輸入量的小樣本數據后，通過最大熵原理計算出輸入量的概率密度函數，可以看出輸入量的概率分布不一定為對稱分布，利用舍選法抽樣產生M個偽隨機數，解決了輸入量小樣本的問題。

將在輸入量概率密度函數的約束下產生的偽隨機數代入測量模型，得到推力矢量參數的概率分布，輸出量的概率分布不一定是正態分布或t分布，通過M個輸出量的值求出不確定度以及給定包含概率下的包含區間。

3.3.1 MCM的適用條件

基于MCM的測量不確定度評估相比GUM法具有以下優點[20]：

1)對模型沒有非線性的限制；

2)不受輸入量的相關性和模型復雜性的影響；

3)不受輸入量分布的影響；

4)不用假設被測量的分布；

5)不用計算偏導數和有效自由度。

3.4 GUM法與MCM的對比驗證

MCM提供了驗證GUM法的方法，當兩者一致時，GUM法仍然作為測量不確定度的主要方法，當兩者不一致時，應采用MCM。

驗證方法：

1)應用GUM法得到輸出量的概率p的包含區間y±Up；

2)通過運用自適應蒙特卡洛法獲得輸出量的標準不確定度u(y)以及概率對稱或最短包含區間的端點值ylow和yhigh；

3)確定由GUM法及MCM獲得的包含區間在約定的數值容差dt下是否一致，確定兩個包含區間的各自端點的絕對偏差：

dlow=|y-Up-ylow|

(19)

dhigh=|y+Up-yhigh|

(20)

如果dlow≤dt，dhigh≤dt，則GUM法通過驗證。

以α，δ為例：

1)推力斜偏角α。有效數字保留3位，數值容差dt=0.5×10-5，dlow=0.042 3-0.039 6=0.000 27>dt，dhigh=0.024 7-0.026 5=0.001 8>dt，GUM法沒有通過MCM的驗證。

2)推力偏移δ。有效數字保留3位，數值容差dt=0.5×10-3，dlow=0.968-0.785=0.183>dt，dhigh=0.672-0.492=0.18>dt，GUM法沒有通過MCM的驗證，相比MCM法較為保守。

4 結束語

對于壓電式推力矢量不確定度評估，數學模型使用了復雜的非線性模型，輸入量的概率分布不是對稱分布，使用GUM法存在一定的局限性，MCM作為GUM法的補充，有更廣的適用范圍，依據各個輸入量的概率密度函數進行隨機抽樣，計算得出M個結果，不僅解決了小樣本測量數據少的問題，減少了大量人力物力的投入，還可以得出不用近似的包含區間，可以驗證GUM法的有效性，可以作為發動機試驗中不確定度評估的有效方法。