基于線性非線性回歸分析的防空導(dǎo)彈費(fèi)用建模*

趙曰強(qiáng)，安實(shí)，麥強(qiáng)，于道林，張耀東

(1.哈爾濱工業(yè)大學(xué) 管理學(xué)院，黑龍江 哈爾濱 150001；2.中國航天科工防御技術(shù)研究院，北京 100854；3.北京電子工程總體研究所，北京 100854)

0 引言

研究防空導(dǎo)彈武器系統(tǒng)費(fèi)用時(shí)，導(dǎo)彈的生產(chǎn)費(fèi)用單價(jià)最受關(guān)注。這是因?yàn)樵谖淦飨到y(tǒng)中導(dǎo)彈的技術(shù)性能指標(biāo)可比性強(qiáng)，而地面系統(tǒng)由于使用要求不同，配置繁簡差別顯著，費(fèi)用的統(tǒng)計(jì)測(cè)算離散很大；從整個(gè)壽命周期費(fèi)用來看，生產(chǎn)費(fèi)用最穩(wěn)定，研制階段由于受到技術(shù)基礎(chǔ)、研制周期、試驗(yàn)鑒定次數(shù)的影響費(fèi)用離散很大；另外，現(xiàn)實(shí)中很多防空導(dǎo)彈武器系統(tǒng)的研制生產(chǎn)，一般將導(dǎo)彈和雷達(dá)制導(dǎo)等地面系統(tǒng)由不同的承包商分開承擔(dān)，導(dǎo)彈相對(duì)獨(dú)立性強(qiáng)。綜上所述，從導(dǎo)彈的生產(chǎn)費(fèi)用單價(jià)入手易于得到可靠的費(fèi)用模型。

1 費(fèi)用變量的分析與選擇

費(fèi)用變量是指影響導(dǎo)彈費(fèi)用的敏感量[1]。防空導(dǎo)彈武器系統(tǒng)性能指標(biāo)較多，各指標(biāo)對(duì)費(fèi)用都有一定的影響，相互間還存在著復(fù)雜的關(guān)聯(lián)性[2]。下面對(duì)主要技術(shù)性能參數(shù)進(jìn)行相關(guān)性分析與選擇。

首先分析導(dǎo)彈殺傷空域?qū)M(fèi)用的影響。一般來說在有效載荷(一般指戰(zhàn)斗部)給定的情況下，射程增大往往造成發(fā)射質(zhì)量增大，進(jìn)而引起費(fèi)用的增加。另外，在射程一定的情況下，提高推進(jìn)劑的性能雖然可以降低導(dǎo)彈發(fā)射質(zhì)量，但推進(jìn)劑的費(fèi)用必定隨著增大。提高制導(dǎo)精度可以降低有效載荷，進(jìn)而大幅降低發(fā)射質(zhì)量，但制導(dǎo)系統(tǒng)的成本會(huì)大幅度提升，這在一定程度上會(huì)抵消導(dǎo)彈發(fā)射質(zhì)量降低帶來的好處。可以看出，在一定制導(dǎo)體制下，射程是影響費(fèi)用的主要因素。

然后分析導(dǎo)彈質(zhì)量對(duì)費(fèi)用的影響。發(fā)射質(zhì)量決定了導(dǎo)彈的大小，從而也決定了所需的動(dòng)力消耗。發(fā)射質(zhì)量越大，費(fèi)用越大；戰(zhàn)斗部質(zhì)量決定了導(dǎo)彈的殺傷威力，在炸藥類型相同的情況下，威力增大，戰(zhàn)斗部質(zhì)量增大，對(duì)制導(dǎo)精度的要求降低?？梢钥闯?，戰(zhàn)斗部質(zhì)量的變化將從導(dǎo)彈質(zhì)量和制導(dǎo)精度2個(gè)方面影響導(dǎo)彈的費(fèi)用，其對(duì)經(jīng)費(fèi)的影響不是單一的增減關(guān)系。

此外分析導(dǎo)彈飛行速度對(duì)費(fèi)用的影響。導(dǎo)彈飛行速度增大，可縮短射擊周期，提高火力密度，增大有效發(fā)射區(qū)和殺傷區(qū)。同時(shí)，導(dǎo)彈飛行速度的提高對(duì)導(dǎo)彈控制能力、結(jié)構(gòu)、防熱水平等都提出了更高的要求，進(jìn)而導(dǎo)致導(dǎo)彈費(fèi)用的增加。

導(dǎo)彈殺傷概率直接取決于導(dǎo)彈的制導(dǎo)精度、引戰(zhàn)配合和殺傷威力。其中，制導(dǎo)精度對(duì)導(dǎo)彈殺傷概率的影響至關(guān)重要。

導(dǎo)彈的可靠性和維修性相關(guān)費(fèi)用可以按導(dǎo)彈武器系統(tǒng)總費(fèi)用的一定比例來概算[3]。

通過以上分析，結(jié)合對(duì)導(dǎo)彈武器系統(tǒng)費(fèi)用敏感的幾個(gè)主要參數(shù)(殺傷空域、導(dǎo)彈質(zhì)量、飛行速度、戰(zhàn)斗部質(zhì)量、制導(dǎo)精度、目標(biāo)通道等)，選取以下參數(shù)作為研究影響導(dǎo)彈單價(jià)的自變量：最大射程、發(fā)射質(zhì)量、最大飛行速度、戰(zhàn)斗部質(zhì)量以及制導(dǎo)誤差。

在較早的模型中，自變量多選用發(fā)射質(zhì)量、發(fā)動(dòng)機(jī)裝藥、戰(zhàn)斗部質(zhì)量、最大射程等組合[4-5]。這是由于20世紀(jì)80年代以前的防空導(dǎo)彈多是第1代、第2代的產(chǎn)品，技術(shù)水平接近，例如：均采用指令制導(dǎo)、脈沖多普勒雷達(dá)。后來隨著技術(shù)的發(fā)展，出現(xiàn)了相控陣?yán)走_(dá)，防空導(dǎo)彈可以打擊多目標(biāo)。于是出現(xiàn)了以最大射程和多目標(biāo)數(shù)作為自變量的組合。但是現(xiàn)今的導(dǎo)彈制導(dǎo)控制技術(shù)得到了大幅度的提升，TVM制導(dǎo)、主動(dòng)雷達(dá)末制導(dǎo)、光學(xué)制導(dǎo)和使用直接力控制使得制導(dǎo)精度越來越高，甚至具備了直接碰撞的能力。所以傳統(tǒng)的變量參數(shù)(如戰(zhàn)斗部質(zhì)量mz)已經(jīng)難以包含影響導(dǎo)彈費(fèi)用的全部主要因素；因此，本文引入了新的變量制導(dǎo)誤差ρ。另外，多目標(biāo)能力在全系統(tǒng)費(fèi)用的估算中將是關(guān)鍵因素之一，系統(tǒng)費(fèi)用建模時(shí)要予以關(guān)注[6]。

2 導(dǎo)彈生產(chǎn)費(fèi)用單價(jià)的線性建模

關(guān)于防空導(dǎo)彈的生產(chǎn)費(fèi)用單價(jià)和性能參數(shù)，數(shù)據(jù)來源較多，數(shù)據(jù)的可信度較高，因此可采用參數(shù)法。具體而言，首先從已經(jīng)選出的有代表性的性能參數(shù)中確定對(duì)導(dǎo)彈費(fèi)用敏感的自變量，例如：最大射程R，發(fā)射質(zhì)量m，導(dǎo)彈最大飛行速度vm以及導(dǎo)彈制導(dǎo)誤差ρ。然后，收集整理類似其他項(xiàng)目的歷史數(shù)據(jù)。最后利用統(tǒng)計(jì)回歸分析技術(shù)建立起費(fèi)用與自變量間的關(guān)系。

假定導(dǎo)彈生產(chǎn)費(fèi)用單價(jià)y(Cmp)與各參數(shù)x1～4(R，m，vm，ρ)之間具有如下的線性函數(shù)關(guān)系：

y=b0+b1x1+b2x2+b3x3+b4x4+e，

(1)

式中：e為服從獨(dú)立正態(tài)的隨機(jī)變量。

運(yùn)用多元線性回歸分析來確定式(1)中的待定參數(shù)b0，b1，b2，b3，b4。

可見對(duì)x1～4(R，m，vm，ρ)的n組不同取值(R1，m0，vm1，ρ1)，(R2，m2，vm2，ρ2)，…，(Rn，mn，vmn，ρn)；有對(duì)應(yīng)的y(Cmp)的值Cmp1，Cmp2，…，Cmpn。

(2)

式中：X為自變量矩陣；Y為因變量；B為回歸系數(shù)；E為隨機(jī)誤差，由式(1)得

Y=XB+BE，

(3)

(4)

變換得

(5)

假定(X′X)的逆矩陣(X′X)-1存在，則得到式(5)回歸問題的最小二乘解為

(6)

2.1 大中型導(dǎo)彈生產(chǎn)費(fèi)用單價(jià)線性模型的建立

建模用到的樣本數(shù)據(jù)真實(shí)有效性是取得良好統(tǒng)計(jì)結(jié)果的關(guān)鍵，參照文獻(xiàn)[7-11]等，查找分析幾個(gè)國內(nèi)外最新防空導(dǎo)彈型號(hào)的性能參數(shù)和費(fèi)用數(shù)據(jù)，并參照美國國防部公布的陸軍緊縮指數(shù)，對(duì)于不同時(shí)間價(jià)格水平進(jìn)行比較，整理列于表1中。

由表1中導(dǎo)彈主要性能技術(shù)參數(shù)作為自變量x1～4(R，m，vm，ρ)，導(dǎo)彈生產(chǎn)費(fèi)用單價(jià)y(Cmp)作為因變量，進(jìn)行多元線性回歸分析，運(yùn)用Matlab運(yùn)算得：

表1 幾種防空導(dǎo)彈費(fèi)用數(shù)據(jù)

注：Cmp為折算的2010年數(shù)據(jù)，單位人民幣萬元

Y= (45.0，49.5，55.6，930.0，157.5，190.0，229.3，

192.5，264.0，316.0，250.0，430.0，690.0，885.0，

856.0，865.0，1 035.0，1 620.0，1 500.0)T.

(7)

將式(7)代入式(2)，得

(8)

代入式(1)，則得回歸費(fèi)用模型：

y= 61.940 2+6.362 1x1-0.126 8x2+

52.780 6x3-13.426 3x4.

(9)

首先利用表1樣本值和計(jì)算值對(duì)模型進(jìn)行驗(yàn)證，由式(9)得出計(jì)算值和歷史數(shù)據(jù)值進(jìn)行比較可以看出，小型導(dǎo)彈的費(fèi)用估算偏差較大，最大到了279.02%，該模型已經(jīng)不能適用；該模型對(duì)大中型導(dǎo)彈的適用性略好。

這里，解決模型對(duì)小型導(dǎo)彈和中大型導(dǎo)彈擬合結(jié)果偏離的問題，美國陸軍導(dǎo)彈系統(tǒng)全壽命費(fèi)用模型中對(duì)于研制費(fèi)用提供了3個(gè)CER解析式，其中第一個(gè)只能適用于便攜導(dǎo)彈的費(fèi)用估計(jì)[12-13]。這也得到一些提示，可以對(duì)小型和大中型導(dǎo)彈的費(fèi)用估計(jì)進(jìn)行分段處理。通過對(duì)導(dǎo)彈數(shù)據(jù)的分析，選用導(dǎo)彈發(fā)射質(zhì)量30 kg作為分段點(diǎn)。

先來看發(fā)射質(zhì)量大于30 kg的導(dǎo)彈費(fèi)用的估計(jì)。將表1中序號(hào)5～19導(dǎo)彈發(fā)射質(zhì)量大于30 kg的導(dǎo)彈主要性能技術(shù)參數(shù)作為自變量，進(jìn)行多元線性回歸，則得大中型導(dǎo)彈的回歸費(fèi)用模型：

y= 192.866 9+5.483 0x1+0.369 2x2+

61.753 7x3-27.243 9x4，m>30 kg.

(10)

同樣先利用樣本值和計(jì)算值對(duì)模型進(jìn)行驗(yàn)證，由式(10)得出計(jì)算值，和樣本值對(duì)比看出，對(duì)于大中型導(dǎo)彈通過多元線性回歸得到的導(dǎo)彈生產(chǎn)單價(jià)的費(fèi)用估算模型并由此而得到的計(jì)算值和歷史數(shù)據(jù)(樣本值)吻合的較好，有一半的數(shù)據(jù)偏差在5%以內(nèi)，最大的偏差僅為21%。

我們進(jìn)一步對(duì)模型的多元回歸進(jìn)行方差分析，分析結(jié)果見表2。

表2 多元線性回歸的方差分析Table 2 Variance analysis of multivariatelinear regression

相關(guān)系數(shù)值為

R2=3 167 900/3 261 700=0.971 3.

(11)

調(diào)整相關(guān)系數(shù)值為

Ra2=1-9 376.3/(3 261 700/14)=0.959 8.

(12)

可見有超過95%的費(fèi)用發(fā)生是由選定的4個(gè)變量(R，m，vm，ρ)決定的。

Ra2和R2數(shù)值差別不大，沒有出現(xiàn)Ra2數(shù)值遠(yuǎn)小于R2的情況，說明參數(shù)選擇數(shù)量合理。

查F分布表可得F4，10(0.01)=5.99<84.465 5(統(tǒng)計(jì)量F的值)，同時(shí)檢驗(yàn)的概率P<0.01，所以，假設(shè)H0：b1=b2=b3=b4=0不成立，回歸方程通過了顯著性檢驗(yàn)，即認(rèn)為所得的線性回歸方程有意義。

2.2 小型導(dǎo)彈生產(chǎn)費(fèi)用單價(jià)線性模型的建立

現(xiàn)在將幾種發(fā)射質(zhì)量小于30 kg的導(dǎo)彈生產(chǎn)單價(jià)和主要性能技術(shù)參數(shù)整理見表3。

從對(duì)表2中的數(shù)據(jù)分析可以看出，小型便攜導(dǎo)彈制導(dǎo)體制相近，制導(dǎo)精度相當(dāng)，導(dǎo)彈生產(chǎn)單價(jià)和作戰(zhàn)空域、速度以及戰(zhàn)斗部質(zhì)量關(guān)聯(lián)不大，而和導(dǎo)彈的發(fā)射質(zhì)量關(guān)聯(lián)顯著，通過各種參數(shù)組合為自變量的回歸及方差分析也可以說明這一點(diǎn)。

將表3中導(dǎo)彈發(fā)射質(zhì)量作為自變量，導(dǎo)彈生產(chǎn)單價(jià)作為因變量，進(jìn)行線性回歸，得小型導(dǎo)彈的回歸費(fèi)用模型：

y=-7.489 0+5.564 7x2，m≤30 kg.

(13)

同樣先利用樣本值和計(jì)算值對(duì)模型進(jìn)行驗(yàn)證，由式(13)得出計(jì)算值，和樣本值對(duì)比看出，由導(dǎo)彈生產(chǎn)單價(jià)的費(fèi)用估算模型得到的計(jì)算值和歷史數(shù)據(jù)(樣本值)最大偏離13%，偏差不大。

進(jìn)一步對(duì)回歸模型進(jìn)行方差分析，見表4。

對(duì)α=0.05，查F分布表可得F1，6(0.05)=5.99。

因?yàn)镕>F1，6(α)，同時(shí)，p=0.026 2<0.05，線性回歸方程反映了導(dǎo)彈生產(chǎn)單價(jià)和導(dǎo)彈發(fā)射質(zhì)量的相關(guān)關(guān)系。

2.3 導(dǎo)彈生產(chǎn)費(fèi)用單價(jià)多元線性回歸模型

根據(jù)上面的計(jì)算和仿真，以導(dǎo)彈發(fā)射質(zhì)量30 kg為分界線建立分段函數(shù)，建立起導(dǎo)彈的生產(chǎn)費(fèi)用單價(jià)y(Cmp)的線性模型。

(14)

式中：y(Cmp)為導(dǎo)彈的生產(chǎn)費(fèi)用單價(jià)，萬元(人民幣)；x1(R)為導(dǎo)彈的發(fā)射質(zhì)量，kg；x2(m)為導(dǎo)彈的最大射程，km；x3(vm)為導(dǎo)彈飛行的最高M(jìn)a；x4(ρ)為導(dǎo)彈制導(dǎo)誤差，m。

表3 幾種小型防空導(dǎo)彈費(fèi)用數(shù)據(jù)Table 3 Cost data of several small air defense missiles

注：表中Cmp為折算的2010年數(shù)據(jù)，單位人民幣萬元。

表4 小型導(dǎo)彈多元線性回歸的方差分析Table 4 Variance analysis of multivariate linearregression for small missile

3 導(dǎo)彈生產(chǎn)費(fèi)用單價(jià)的非線性建模

上節(jié)建立起了導(dǎo)彈的生產(chǎn)費(fèi)用單價(jià)與特性參數(shù)之間的線性函數(shù)關(guān)系。采用了對(duì)小型和大中型導(dǎo)彈的費(fèi)用估計(jì)進(jìn)行分段處理，取得了理想的結(jié)果。下面進(jìn)一步開展非線性分析，建立非線性回歸模型，來探求比較形成最佳的解決方案[14]。

基本原理是非線性回歸最小二乘擬合[15]。利用Matlab非線性曲線擬合命令nlinfit函數(shù)，進(jìn)行非線性函數(shù)的系數(shù)用最小二乘法估計(jì)，該函數(shù)使用高斯—牛頓算法，調(diào)用格式為

beta=nlinfit(X,y,@myfun,beta0).

(15)

式中：y為因變量矢量；X為自變量組成的矩陣；beta0為包含系數(shù)初始值的矢量；@myfun參數(shù)為一函數(shù)，具有以下的形式

yhat=myfun(beta,X)，

(16)

式中：beta為系數(shù)矢量；X為自變量矩陣；@myfun為參數(shù)返回一個(gè)擬合y值的yhat矢量。

3.1 建立二次函數(shù)模型

設(shè)定y(Cmp)與各參數(shù)x1～5(R，m，mz，vm，ρ)之間具有如下的函數(shù)關(guān)系：

y= beta1x1+beta2x2+beta3x3+beta4x4+beta5x5+

beta6x12+beta7x1x2+beta8x1x3+beta9x1x4+

beta10x1x5+beta11x22+beta12x2x3+beta13x2x4+

beta14x2x5+beta15x32+beta16x3x4+beta17x3x5+

beta18x42+beta19x4x5+beta20x52.

(17)

這里構(gòu)建的是含有平方項(xiàng)和交互項(xiàng)的二次函數(shù)表達(dá)式。

調(diào)用Matlab中nlinfit函數(shù)進(jìn)行估算，得到beta為系數(shù)矢量，beta(20)為極小值，且超出置信區(qū)間略去，即得到最終的多元二次非線性回歸的解析式：

y= -221.64x1-49.448x2+609.5x3+168.38x4-

144.07x5-0.1807 5x12-0.249 01x1x2+

0.816 3x1x3+72.346x1x4+3.574 5x1x5-

0.035 167x22-0.269 65x2x3+11.512x2x4+

4.598 9x2x5+10.215x32-192.42x3x4-

60.829x3x5-105.38x42+202.59x4x5，

(18)

式中：y(Cmp)為導(dǎo)彈的生產(chǎn)費(fèi)用單價(jià)，萬元(人民幣)；x1(R)為導(dǎo)彈的發(fā)射質(zhì)量，kg；x2(m)為導(dǎo)彈的最大射程，km；x3(mz)為戰(zhàn)斗部的質(zhì)量，kg；x4(vm)為導(dǎo)彈飛行的最高M(jìn)a；x5(ρ)為導(dǎo)彈制導(dǎo)誤差，m。

對(duì)該模型進(jìn)行回歸檢驗(yàn)，由命令函數(shù)[b,bint,r,rint,stats]=regress(y_result,x)，得檢驗(yàn)回歸模型的統(tǒng)計(jì)量：

stats=(0.959 47，86.075，1.074 4e-9，12 430).

(19)

可見，相關(guān)系數(shù)平方值R2=0.9594 7說明回歸方程顯著模型擬合程度較高，F(xiàn)=86.075不超臨界，顯著性概率P=1.074 4e-9<0.05拒零假設(shè)，模型有意義，均方根誤差RMSE=12 430。

采用式(18)，利用表1中的導(dǎo)彈費(fèi)用數(shù)據(jù)擬合結(jié)果為

Y_result= (45，49.5，55.6，93，157.5，190，229.3，

192.5，264，316，250，430，690，885，

856，865，1 035，1 620，1 500).

(20)

將擬合值與樣本值進(jìn)行對(duì)比見圖1。

圖1 二次函數(shù)模型擬合值與樣本值對(duì)比圖Fig.1 Contrast diagram of fitting value and samplevalue of quadratic function model

本模型擬合精度高，用一個(gè)解析式較一致地表達(dá)了因變量和各參數(shù)的關(guān)系，而不需要采用分段表示。相對(duì)線性回歸模型，該模型具有更高的計(jì)算復(fù)雜性，適合使用計(jì)算機(jī)進(jìn)行計(jì)算。

3.2 建立非線性函數(shù)模型

為了進(jìn)一步簡化允許特定一點(diǎn)出現(xiàn)較大偏差(可視孤立樣本予以剔除)的條件下的擬合的函數(shù)。先看下每個(gè)參數(shù)對(duì)性能的影響，確定交互項(xiàng)和單一項(xiàng)，然后遍歷所有指數(shù)的可能，獲得和原始數(shù)據(jù)的最小誤差。建立生產(chǎn)費(fèi)用單價(jià)的模型：

y= -71.589+8.692x1-727.39x4- 9.710 4x10.5x30.5+749.1x20.05x40.9，

(21)

式中：y(Cmp)為導(dǎo)彈的生產(chǎn)費(fèi)用單價(jià)，萬元(人民幣)；x1(R)為導(dǎo)彈的發(fā)射質(zhì)量，kg；x2(m)為導(dǎo)彈的最大射程，km；x3(mz)為戰(zhàn)斗部的質(zhì)量，kg；x4(vm)為導(dǎo)彈飛行的最高M(jìn)a。

采用上式，利用表1中的導(dǎo)彈費(fèi)用數(shù)據(jù)擬合結(jié)果為

Y_result= (61.364，49.017，49.017，73.463，

183.54，193.23，189.77，206.49，228.13，

308.03，281.18，540.08，681.78，911.92，

895 718.29，1 040.4，1 613.6，1 499.2).

(22)

將擬合值與樣本值進(jìn)行對(duì)比如圖2所示。

圖2 非線性模型擬合值與樣本值對(duì)比圖Fig.2 Contrast diagram of fitting value andsample value of nonlinear model

對(duì)該模型進(jìn)行回歸檢驗(yàn)，得檢驗(yàn)回歸模型的統(tǒng)計(jì)量

stats=(0.968 53，110.25，2.032 3e-10，9 557.2).

(23)

可見，相關(guān)系數(shù)平方值R2=0.968 5說明方程回歸顯著，顯著性概率P=2.032 3 e-10<0.05存在非零系數(shù)，說明擬合模型有效。

該模型式(21)和原數(shù)據(jù)總的誤差為540，平均每個(gè)點(diǎn)的誤差為28.4，最大點(diǎn)的誤差為150.5。

本模型相對(duì)二次函數(shù)模型式(18)大為簡化，檢驗(yàn)有效，但是從圖2可知，整體擬合較好，個(gè)別點(diǎn)有一定的偏差(17%以內(nèi))。

4 結(jié)束語

采用線性回歸的方法建立導(dǎo)彈生產(chǎn)費(fèi)用模型式(14)，由于選取的參數(shù)典型且與費(fèi)用相關(guān)大，建立的模型相對(duì)簡單且精度高，易于計(jì)算，且物理意義明確，不足之處是要根據(jù)導(dǎo)彈的總質(zhì)量大小分兩段進(jìn)行擬合。非線性回歸二次函數(shù)建模、非線性建模都取得了很好的結(jié)果，避免了線性回歸模型分段的不足；但計(jì)算略復(fù)雜，適合計(jì)算機(jī)運(yùn)算。二次函數(shù)模型式(18)精度除質(zhì)量小、近程的個(gè)別裝備，整體精度較高，但結(jié)構(gòu)復(fù)雜，幾乎包含了各參數(shù)的平方項(xiàng)和交互項(xiàng)；非線性模型式(21)結(jié)構(gòu)簡單，精度較好，能用一個(gè)式子簡明準(zhǔn)確地表達(dá)了所有防空導(dǎo)彈的生產(chǎn)費(fèi)用單價(jià)。本文建模選取的費(fèi)用變量雖然相關(guān)性較高，但在應(yīng)用中還應(yīng)根據(jù)工程實(shí)際，考慮其他費(fèi)用變量的作用。