鄭 科, 耿衛(wèi)國,朱子環(huán)

(北京航天試驗(yàn)技術(shù)研究所，北京 100074)

0 引言

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

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

1 壓電式推力矢量測量

火箭發(fā)動機(jī)壓電式推力矢量系統(tǒng)的數(shù)學(xué)測量模型如圖1所示[7]。

圖1 某壓電式推力矢量測量系統(tǒng)數(shù)學(xué)模型

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

F為姿軌控發(fā)動機(jī)點(diǎn)火產(chǎn)生的空間矢量推力；

Fx，F(xiàn)y，F(xiàn)z分別為矢量推力F在作用點(diǎn)(δy,δz)的三向分力；

Fx1，F(xiàn)x2，F(xiàn)x3，F(xiàn)x4分別為通過測力平臺上的4個排列成正方形分布的三維力傳感器實(shí)際測得的X方向的4個分力。同樣，F(xiàn)y1，F(xiàn)y2，F(xiàn)y3，F(xiàn)y4和Fz1，F(xiàn)z2，F(xiàn)z3，F(xiàn)z4分別為通過傳感器實(shí)際測得的Y方向和Z方向分力；

Mo為火箭發(fā)動機(jī)矢量推力F對測量平臺坐標(biāo)系中心產(chǎn)生的總力距；

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

a為發(fā)動機(jī)測力平臺上三維力傳感器與坐標(biāo)軸之間的距離；

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

建立火箭發(fā)動機(jī)推力矢量參數(shù)的測量模型：

發(fā)動機(jī)的三項推力為:

(1)

推力矢量力矩為:

(2)

則得:

推力斜偏角α:

(3)

側(cè)向力方位角:

(4)

推力偏移:

(5)

推力偏移方位角:

(6)

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

2.1 蒙特卡洛法原理

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

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

(7)

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

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

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

MCM評估過程：

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

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

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

4)從各個輸入量Xi的概率密度函數(shù)中隨機(jī)抽樣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的期望值以及標(biāo)準(zhǔn)差，求出估計值y和不確定度μ(y)；

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

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

圖3 蒙特卡洛法流程圖

2.2 蒙特卡洛法的實(shí)現(xiàn)

2.2.1 概率密度函數(shù)的確定

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

2.2.2 抽樣仿真次數(shù)M

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

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

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

(8)

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

(9)

根據(jù)自由度的計算公式可得:

(10)

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

根據(jù)上述兩種方法得到的計算結(jié)果，結(jié)合不確定度評估的實(shí)際經(jīng)驗(yàn)，一般可以取抽樣仿真次數(shù)M為106。

2.2.3 舍選法抽樣

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

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

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

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

(11)

2.2.4 蒙特卡洛法仿真結(jié)果

(12)

(13)

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

3.1 測量數(shù)據(jù)

根據(jù)某次推力試車數(shù)據(jù)推導(dǎo)得出推力矢量參數(shù)如表1所示。

表1 某次試車推力矢量參數(shù)

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

表2 壓電式推力矢量系統(tǒng)X向靜態(tài)標(biāo)定

3.2 GUM法不確定的評估

計算參數(shù)合成不確定度一般表達(dá)式為:

(14)

式中,f為被測量y與直接測得量xi的函數(shù)關(guān)系。

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

根據(jù)合成標(biāo)準(zhǔn)不確定度推導(dǎo)：

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

(15)

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

(16)

由前面章節(jié)推導(dǎo)公式得出推力矢量參數(shù)的合成不確定度。例如:

1)α的合成不確定度:

(17)

2)δ的合成不確定度:

(18)

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

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

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

3.2.1 GUM法的適用條件

在實(shí)際應(yīng)用中，GUM法適用于以下條件：

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

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

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

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

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

3.3 MCM的軟件實(shí)現(xiàn)

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

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

軟件具體實(shí)施方式：

1)將發(fā)動機(jī)試驗(yàn)數(shù)據(jù)導(dǎo)入軟件界面；

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

3)確定每個輸入量的數(shù)據(jù)邊界；

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

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

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

將表3X方向的靜態(tài)標(biāo)定線性差值，以及表4Y方向、表5Z方向的靜態(tài)標(biāo)定線性差值導(dǎo)入不確定度評估軟件中，如圖4所示，不確定度評估結(jié)果如表6所示。

表3 壓電式推力矢量系統(tǒng)Y向靜態(tài)標(biāo)定

表4 壓電式推力矢量系統(tǒng)Z向靜態(tài)標(biāo)定

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

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

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

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

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

3.3.1 MCM的適用條件

基于MCM的測量不確定度評估相比GUM法具有以下優(yōu)點(diǎn)[20]：

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

2)不受輸入量的相關(guān)性和模型復(fù)雜性的影響；

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

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

5)不用計算偏導(dǎo)數(shù)和有效自由度。

3.4 GUM法與MCM的對比驗(yàn)證

MCM提供了驗(yàn)證GUM法的方法，當(dāng)兩者一致時，GUM法仍然作為測量不確定度的主要方法，當(dāng)兩者不一致時，應(yīng)采用MCM。

驗(yàn)證方法：

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

2)通過運(yùn)用自適應(yīng)蒙特卡洛法獲得輸出量的標(biāo)準(zhǔn)不確定度u(y)以及概率對稱或最短包含區(qū)間的端點(diǎn)值ylow和yhigh；

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

dlow=|y-Up-ylow|

(19)

dhigh=|y+Up-yhigh|

(20)

如果dlow≤dt，dhigh≤dt，則GUM法通過驗(yàn)證。

以α，δ為例：

1)推力斜偏角α。有效數(shù)字保留3位，數(shù)值容差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的驗(yàn)證。

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

4 結(jié)束語

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