小樣本下煤油加注溫升模型參數(shù)估計(jì)方法

劉巾杰，許琪琪，扶山川，竇天恒

(西昌衛(wèi)星發(fā)射中心，海南 文昌 571300)

0 引言

煤油作為目前液體火箭的主要推進(jìn)劑之一，其溫度波動(dòng)會(huì)引起物性參數(shù)變化，對(duì)火箭發(fā)動(dòng)機(jī)的性能產(chǎn)生影響。煤油加注系統(tǒng)是大推力運(yùn)載火箭動(dòng)力系統(tǒng)的重要組成部分，其主要功能是連續(xù)、可靠、準(zhǔn)確地向箭上貯箱加注煤油[1]。在加注過(guò)程中，受到環(huán)境溫度、加注管路、加注泵工作狀態(tài)等因素的影響，煤油溫度會(huì)發(fā)生顯著變化[2-3]。因此，需要對(duì)煤油初始溫度進(jìn)行調(diào)整，控制煤油進(jìn)箭溫度。

采用準(zhǔn)確的溫升參數(shù)，建立科學(xué)的溫升模型，才能對(duì)煤油溫升過(guò)程進(jìn)行精確模擬，進(jìn)而根據(jù)最終溫度要求反算出最優(yōu)的煤油初始溫度[2]。目前靶場(chǎng)采用的溫升參數(shù)，皆為確定的經(jīng)驗(yàn)參數(shù)。然而以往試驗(yàn)任務(wù)與新型運(yùn)載火箭靶場(chǎng)實(shí)際加注設(shè)備、停放環(huán)境等差別較大，若仍然使用經(jīng)驗(yàn)參數(shù)進(jìn)行計(jì)算，會(huì)致使計(jì)算溫度與煤油經(jīng)過(guò)加注溫升后的實(shí)際溫度之間存在不可忽略的相對(duì)誤差。因此有必要根據(jù)靶場(chǎng)實(shí)測(cè)溫度數(shù)據(jù)，對(duì)煤油加注溫升參數(shù)進(jìn)行估計(jì)。目前，我國(guó)靶場(chǎng)進(jìn)行液體火箭飛行試驗(yàn)的次數(shù)極少，如何利用僅有的幾次實(shí)測(cè)數(shù)據(jù)樣本，設(shè)計(jì)更為有效的估計(jì)方法，對(duì)溫升參數(shù)進(jìn)行估計(jì)，使煤油溫度更好地滿足發(fā)動(dòng)機(jī)工作要求，對(duì)提升火箭發(fā)動(dòng)機(jī)整體性能有著重要意義。

1 加注溫升模型

常溫煤油主要采用泵壓式加注[1]。加注溫升是指加注結(jié)束時(shí)火箭貯箱中的煤油溫度(不考慮溫度分層)與加注前庫(kù)區(qū)儲(chǔ)罐中煤油溫度之差，其計(jì)算模型為

Tj=Tk+f1(Ta-Tk)+f0+ε

(1)

式中：Tj為煤油加注結(jié)束時(shí)溫度；Tk為煤油庫(kù)房?jī)?chǔ)罐內(nèi)溫度；f1為加注溫變系數(shù)；Ta為煤油加注開始時(shí)刻、結(jié)束時(shí)刻環(huán)境溫度的平均值；f0為加注泵溫升；ε為測(cè)量誤差。

在實(shí)際煤油加注前，根據(jù)目標(biāo)溫度，依照給定的加注模型，對(duì)初始加注溫度進(jìn)行反算。而正式加注時(shí)，是給定初始加注溫度的情況下，煤油從儲(chǔ)罐加注到火箭貯箱的過(guò)程中，由于加注泵加壓、環(huán)境熱交換等多種原因，溫度自然升高，獲得加注最終溫度。事實(shí)上，加注溫升模型，就是對(duì)煤油在加注過(guò)程中溫度升高過(guò)程的一種模擬。對(duì)同一個(gè)溫升模型，從同一個(gè)初始溫度出發(fā)，其溫升參數(shù)越精確，解算出來(lái)的加注結(jié)束溫度就與實(shí)際的加注結(jié)束溫度最接近。圖1為A,B兩次任務(wù)加注溫度數(shù)據(jù)。其中，A任務(wù)數(shù)據(jù)最大相對(duì)誤差達(dá)到1.099 5 ℃，B任務(wù)數(shù)據(jù)最大相對(duì)誤差達(dá)到4.557 9 ℃，均超出了煤油溫度計(jì)算要求精度(要求為±1 ℃)。

圖1 兩次任務(wù)加注溫度數(shù)據(jù)Fig.1 Filling temperature of two tests

不妨假設(shè)式(1)所示的模型是正確而完備的，其中Tj，Tk，Ta可以通過(guò)實(shí)測(cè)獲得，可見(jiàn)造成計(jì)算結(jié)果相對(duì)誤差的根本在于使用的經(jīng)驗(yàn)參數(shù)f1，f0不夠準(zhǔn)確。

2 溫升模型的貝葉斯回歸分析

xi=Ta-Tk

令

Y=[y1,y2,…,ym]T

β=[f1,f0]T

將式(1)進(jìn)一步寫成

Y=Xβ+ε

(2)

由式(2)可見(jiàn)，這是一個(gè)一元回歸模型[4]。

采用最小二乘求解，有[5]

(3)

小樣本條件下，貝葉斯法進(jìn)行回歸分析是一種較好的估計(jì)方法[9-13]。貝葉斯方法通過(guò)將估計(jì)參數(shù)的先驗(yàn)信息加入回歸分析中，能夠有效降低測(cè)量誤差對(duì)參數(shù)估計(jì)的影響，使小樣本條件下仍然可以得到參數(shù)的合理估計(jì)[14-15]。

不妨假設(shè)測(cè)量的隨機(jī)誤差ε是一個(gè)高斯白噪聲，式(2)中有：E(ε)=0，Cov(ε)=σ2Im。并且有先驗(yàn)分布：β～N(E(β),σ2I2)，其中E(β)=μ，(μ為已知超參數(shù))。在煤油加注溫升模型中，μ為給定的經(jīng)驗(yàn)參數(shù)：加注溫變系數(shù)f1和加注泵溫升f0，記為μ=[f1,f0]T。

在式(2)所示的回歸模型中參數(shù)的最小風(fēng)險(xiǎn)貝葉斯線性無(wú)偏估計(jì)(簡(jiǎn)稱MRBLUE)定義如下[8]：

其中

(4)

b=(Im-AX)μ

(5)

計(jì)算貝葉斯風(fēng)險(xiǎn)

E[A(Y-Xμ)-(β-μ)]TD[A(Y-Xμ)-(β-μ)]=

Etr{D[A(Y-Xμ)-(β-μ)][A(Y-Xμ)-(β-μ)]T}=

σ2tr{DA(Im+XXT)AT+DAX-DXTAT}

(P+BCBT)-1=P-1-P-1B(BTP-1B+C-1)-1BTP-1

即

A=(XXT+Im)-1XT

(6)

所以β的最小風(fēng)險(xiǎn)貝葉斯線性無(wú)偏估計(jì)為

(7)

對(duì)上式進(jìn)行進(jìn)一步梳理，有

(8)

(9)

3 加注溫升參數(shù)估計(jì)

圖2 基于A試驗(yàn)任務(wù)數(shù)據(jù)的加注溫度對(duì)比Fig.2 Comparison of filling temperature based on test A

3種方法測(cè)得的參數(shù)中，由貝葉斯回歸參數(shù)得到的加注結(jié)束溫度與實(shí)測(cè)溫度最為相近。進(jìn)一步計(jì)算3種方法得到的加注結(jié)束溫度與實(shí)測(cè)溫度的誤差，如表1所示。

表1 基于A樣本數(shù)據(jù)的誤差對(duì)比

經(jīng)過(guò)貝葉斯回歸估計(jì)參數(shù)后得到的加注結(jié)束溫度與實(shí)測(cè)溫度最為吻合，相對(duì)經(jīng)驗(yàn)參數(shù)的結(jié)果改善了53.89%。最小二乘結(jié)果與經(jīng)驗(yàn)參數(shù)結(jié)果相似，但最小二乘方法受測(cè)量誤差影響，估算得到的加注泵溫升f0=-6.579 9<0，與實(shí)際情況不符，可見(jiàn)小樣本條件下最小二乘方法并不適用。

將從此次加注過(guò)程中計(jì)算得到的貝葉斯回歸參數(shù)用于該型火箭飛行試驗(yàn)任務(wù)B的加注過(guò)程，使用煤油初始加注溫度計(jì)算得到加注結(jié)束溫度，聯(lián)合使用實(shí)驗(yàn)測(cè)定參數(shù)算出的結(jié)果與實(shí)測(cè)溫度進(jìn)行比較，結(jié)果如圖3所示。

圖3 基于B試驗(yàn)任務(wù)數(shù)據(jù)的加注溫度對(duì)比Fig.3 Comparison of filling temperature based on test B

進(jìn)一步計(jì)算三種方法得到的加注結(jié)束溫度與實(shí)測(cè)溫度的誤差，如表2所示。經(jīng)過(guò)貝葉斯估計(jì)參數(shù)后得到的加注結(jié)束溫度與實(shí)測(cè)溫度最為吻合，相對(duì)經(jīng)驗(yàn)參數(shù)的結(jié)果改善了65.12%。相對(duì)于A樣本數(shù)據(jù)，B樣本數(shù)據(jù)體現(xiàn)出來(lái)的實(shí)驗(yàn)測(cè)定結(jié)果誤差更小，是因?yàn)楸疚闹惺褂玫腁任務(wù)給定的經(jīng)驗(yàn)參數(shù)，其本身就是針對(duì)B樣本中擬合不好的問(wèn)題進(jìn)行修正過(guò)的結(jié)果，所以采用此參數(shù)計(jì)算A任務(wù)數(shù)據(jù)效果要優(yōu)于B任務(wù)數(shù)據(jù)。同時(shí)，可以看到，直接根據(jù)靶場(chǎng)小樣本數(shù)據(jù)進(jìn)行貝葉斯回歸修正后的參數(shù)，其擬合結(jié)果在A、B數(shù)據(jù)中皆優(yōu)于經(jīng)驗(yàn)參數(shù)的擬合結(jié)果。

表2 基于B樣本數(shù)據(jù)的誤差對(duì)比

上述結(jié)果說(shuō)明利用歷史數(shù)據(jù)進(jìn)行貝葉斯回歸得到的參數(shù)與實(shí)驗(yàn)測(cè)定參數(shù)相比，能更好地?cái)M合加注過(guò)程的煤油溫升。據(jù)此驗(yàn)證了貝葉斯回歸進(jìn)行參數(shù)估計(jì)，改進(jìn)未來(lái)煤油加注溫度計(jì)算這一方法的可行性和科學(xué)性。

值得注意的是，貝葉斯回歸不僅在小樣本情況下實(shí)現(xiàn)科學(xué)合理的參數(shù)估計(jì)，它同樣可以在樣本量增加的情況下，提高自身的估計(jì)精度。圖4所示為使用A，B兩次樣本數(shù)據(jù)進(jìn)行估計(jì)的參數(shù)，與只使用A一次樣本數(shù)據(jù)進(jìn)行估計(jì)的參數(shù)，對(duì)A加注結(jié)束溫度進(jìn)行擬合解算的效果對(duì)比。

圖4 樣本數(shù)量對(duì)參數(shù)估計(jì)效果的影響Fig.4 The influence of sample size to estimation result

使用兩次數(shù)據(jù)樣本估計(jì)的參數(shù)，解算出的加注結(jié)束溫度與實(shí)際測(cè)量溫度的均方誤差為4.056 6，與一次數(shù)據(jù)樣本估計(jì)的參數(shù)得到的均方誤差5.937 5(見(jiàn)表1)相比，有了進(jìn)一步改善。

4 結(jié)語(yǔ)

為適應(yīng)未來(lái)火箭飛行試驗(yàn)更大推力、更高精度的要求[15]，精準(zhǔn)加注是靶場(chǎng)必須解決的技術(shù)問(wèn)題。在最快的時(shí)間內(nèi)估計(jì)出最精確的參數(shù)，建立最精確的推進(jìn)劑溫升模型，是其中的關(guān)鍵環(huán)節(jié)。利用僅有的幾次飛行試驗(yàn)數(shù)據(jù)，本文采用貝葉斯回歸分析方法，很好地解決了靶場(chǎng)小樣本下煤油加注溫升參數(shù)估計(jì)問(wèn)題。相比于目前使用的經(jīng)驗(yàn)參數(shù)和經(jīng)典的最小二乘方法，貝葉斯回歸能夠綜合經(jīng)驗(yàn)參數(shù)先驗(yàn)信息和靶場(chǎng)實(shí)地實(shí)測(cè)數(shù)據(jù)進(jìn)行科學(xué)估計(jì)，有效抑制測(cè)量誤差的影響，顯著提高估計(jì)精度。隨著液體火箭飛行試驗(yàn)任務(wù)的逐漸增多，可用的數(shù)據(jù)樣本容量相應(yīng)增大，貝葉斯回歸分析的參數(shù)估計(jì)精度會(huì)得到進(jìn)一步提升，為精準(zhǔn)加注提供更為堅(jiān)實(shí)的理論基礎(chǔ)。