相關(guān)性和顯著性結(jié)合的制導(dǎo)工具誤差系數(shù)分離方法

魏宗康，高榮榮，趙 友

（北京航天控制儀器研究所，北京 100854）

慣性制導(dǎo)工具誤差分離精度對于提升彈道導(dǎo)彈命中精度有著舉足輕重的意義。在利用彈道導(dǎo)彈遙外測數(shù)據(jù)分離慣性測量系統(tǒng)誤差系數(shù)時，可采用遙外測速度誤差作為觀測量，其優(yōu)點是速度誤差反映了加速度計組合和陀螺儀組合的測量誤差，另外一個優(yōu)點是建立速度環(huán)境函數(shù)矩陣后可直接通過解方程求解誤差系數(shù)，過程中沒有微分計算。

基于速度環(huán)境函數(shù)的遙外測誤差分離時，首先要確定誤差模型的結(jié)構(gòu)。彈道導(dǎo)彈飛行軌跡的特點是慣性系統(tǒng)只工作于主動段，而不能全姿態(tài)任意方位都存在大過載或大機動，這就決定了選擇的誤差模型結(jié)構(gòu)矩陣中部分系數(shù)之間具有相關(guān)性，而最小二乘法對于強相關(guān)結(jié)構(gòu)矩陣有較低的適應(yīng)能力，會導(dǎo)致分離的誤差系數(shù)偏離真值較大。為此，如何在強相關(guān)條件的約束下，實現(xiàn)基于彈道導(dǎo)彈遙外測數(shù)據(jù)的慣性測量系統(tǒng)誤差系數(shù)的精確分離是一個難題。為解決該問題，文獻[1]采用靈敏度理論對制導(dǎo)工具誤差進行分離，通過比較兩個誤差系數(shù)的靈敏度曲線變化趨勢的一致性來判斷兩個誤差系數(shù)的相關(guān)性，但是此方法的缺點是沒有具體的量化指標(biāo)。文獻[2]和文獻[3]分別提出了主成分估計、偏最小二乘估計、嶺估計法、粒子群算法等不同的估計方法，這些方法的核心思想是在模型不降階的前提下估計出系數(shù)的近似解，但這些方法帶來的問題就是估計有偏差，沒有從根本上解決慣性制導(dǎo)工具誤差系數(shù)的精確分離問題。文獻[4]采用顯著性分析方法對石英加速度計的二次項誤差系數(shù)進行分離，但是并沒有考慮誤差系數(shù)之間的相關(guān)性問題，分離方法缺乏嚴(yán)謹(jǐn)性。

針對該問題本文提出了一種把相關(guān)性檢驗和顯著性檢驗分階段結(jié)合起來的慣性制導(dǎo)工具誤差分離方法，重點不是在于估計結(jié)構(gòu)矩陣相關(guān)時的每個系數(shù)，而是估計不同相關(guān)系數(shù)之間組合后的值，采用該方法后分離的結(jié)果更準(zhǔn)確，也較符合工程實際。

1 慣性測量系統(tǒng)高階誤差模型

加速度計組合的高階誤差模型分別為[5-6]：

式(1)中，k0x、k0y、k0z為x、y、z加速度計的零偏，單位為g；δkx、δky、δkz為x、y、z加速度計的線性度誤差；δKax、δKay、δKaz為x、y、z加速度計的線性度不對稱性誤差；k yx、k zx、k xy、k zy、k xz、kyz為x、y、z加速度計的安裝誤差角，單位為rad；K2x、K2y、K2z為x、y、z加速度計的二次項誤差系數(shù)，單位為g/g2；δK2x、δK2y、δK2z為x、y、z加速度計的奇二次項誤差系數(shù)，單位為g/g2；K xxy、K xxz、K xyz、Kyxy、Kyxz、Kyyz、Kzxy、Kzxz、Kzyz為x、y、z加速度計的交叉耦合項誤差系數(shù)，單位為g/g2；K3x、K3y、K3z分別為x、y、z加速度計的三次項誤差系數(shù)，單位為g/g3。

陀螺儀組合的高階誤差模型為[7-8]

式(2)中，DFx、DFy、DFz為x、y、z陀螺儀的常值漂移，單位為°/h；D1x、D1yD1z、D2x、D2y、D2z、D3x、D3y、D3z為x、y、z陀螺儀的一次項漂移，單位為°/h/g；D4x、D4y、D4z、D5x、D5y、D5z、D6x、D6y、D6z為x、y、z陀螺儀的二次項漂移，單位為°/h/g2；D7x、D7y、D7z、D8x、D8y、D8z、D9x、D9y、D9z為x、y、z陀螺儀的交叉耦合二次項漂移，單位為°/h/g2。

2 遙外測速度誤差

對于在彈道導(dǎo)彈主動段關(guān)機點就結(jié)束工作的慣性系統(tǒng)來說，由于時間較短，可在誤差分析時忽略各反饋回路的影響，簡化的導(dǎo)航誤差傳播流程見圖1。

圖1 基于地球坐標(biāo)系的導(dǎo)航誤差傳播開環(huán)流程圖Fig.1 Open-loop flow chart of navigation error transmission based on earth coordinate system

采用上述簡化誤差模型的優(yōu)點是可直接對姿態(tài)微分方程、速度微分方程和位置微分方程進行積分，對應(yīng)列寫出各項慣性儀表誤差系數(shù)的環(huán)境函數(shù)矩陣，采用最小二乘法直接求解誤差系數(shù)，該求解方法相對簡單[9]。圖1 中，

Aε——由陀螺漂移引起姿態(tài)角速度誤差時的變換矩陣；

Aφ——由姿態(tài)角誤差引起速度誤差變化量時的變化矩陣；

根據(jù)速度誤差傳播方程即可列寫出把陀螺儀、加速度計各項誤差系數(shù)對應(yīng)的速度環(huán)境函數(shù)[10]。速度誤差方程的模型結(jié)構(gòu)為：

上式采用矩陣形式可寫為：

其中，Y表示速度誤差矩陣；C表示速度環(huán)境函數(shù)矩陣；X表示慣性器件的誤差系數(shù)。

3 主成分估計的缺陷

當(dāng)結(jié)構(gòu)矩陣C不為維列滿秩時，則不能用最小二乘法求解[11-12]。所謂主成分估計就是把信息矩陣Φ=CTC中較小的特征值近似簡化為0，然后求解X的過程。

當(dāng)C不為列滿秩時，Φ的特征值λ1、λ2…λm中至少有一個為0。為分析方便，設(shè)λ1、λ2、…、λm為從小到大排序后的結(jié)果，則D=diag（λ1,λ2…λm）可寫為：

其中，ΛB為對角線所有元素非零的對角陣。

根據(jù)矩陣D對應(yīng)的正交變換矩陣表示為：

則有：

由于PA為信息矩陣Φ=CTC的特征值0 對應(yīng)的特征向量，有

對上式可解出：

由此可知，當(dāng)C不為列滿秩時，求解方程Y=CX等價于求解方程：

一種求解X的應(yīng)用較廣但不準(zhǔn)確的方法是，令：

可得到：

因此，采用主成分估計時，仍然需要非相關(guān)處理。只有在對結(jié)構(gòu)矩陣進行降維處理后，才可求得部分系數(shù)的解。主成分辨識法的另外一個缺點是特征值大并不意味著對應(yīng)誤差系數(shù)也顯著，因此，主成分估計方法無法識別誤差系數(shù)在實際使用中的最優(yōu)估計問題[13]。

4 結(jié)合相關(guān)性和顯著性的辨識方法

針對結(jié)構(gòu)矩陣不為列滿秩的情形，本節(jié)的研究思路則是對模型適當(dāng)降階以求出各系數(shù)不同組合的精確解[14]。

4.1 結(jié)構(gòu)矩陣列向量相關(guān)時的系數(shù)計算

在結(jié)構(gòu)矩陣不為列滿秩時，則必然有相應(yīng)的列向量相關(guān)。在求解系數(shù)時，可先找到結(jié)構(gòu)矩陣中最相關(guān)的列后，去掉其中的一列來對該結(jié)構(gòu)矩陣進行簡化[15]。

下面給出相關(guān)性定義，設(shè)有向量Ci（i=1,2…m）為結(jié)構(gòu)矩陣的列矩陣，當(dāng)Ci全不為零時，若有某些非零標(biāo)量αi（i=1,2…m）滿足：

則稱向量C1、C2···Cm線性相關(guān)。

若僅當(dāng)：

時式(15)才成立，則稱向量C1、C2···Cm線性獨立，或不相關(guān)。

為分析方便，設(shè)向量C1、C2···Cj-1、Cj、Cj+1、···Cm線性相關(guān)，且αj非零，則必然存在某些非零標(biāo)量rj,i（i=1,2…j-1，j+1…m）滿足：

把Cj代入線性方程Y=CX中，有：

4.2 辨識具體步驟

顯著性檢驗只是解決了因子與輸出的關(guān)系，但不能解決因子之間的相關(guān)性問題；相關(guān)性檢驗只是解決了因子之間的相關(guān)性問題，而沒有解決因子與輸出的顯著性問題。下面試圖把顯著性檢驗和相關(guān)性檢驗結(jié)合起來以辨識系統(tǒng)的系數(shù)。

在結(jié)構(gòu)矩陣中存在相關(guān)的兩個列時，采用離線最小二乘法時信息矩陣CTC的行列式將近似為0，此時分離的系數(shù)值將嚴(yán)重偏離真實值，甚至無法求解。雖然采用顯著性檢驗可以簡化模型的結(jié)構(gòu)，但這種簡化也存在一個缺點，就是在結(jié)構(gòu)矩陣奇異或近似奇異時有可能把兩個相關(guān)的變量都作為顯著項而保留到模型結(jié)構(gòu)[16]。

為解決奇異線性系統(tǒng)的離線系數(shù)辨識問題，在本節(jié)提出一種綜合相關(guān)性和顯著性的辨識方法。就是對一個考慮到各種因子的線性系統(tǒng)模型結(jié)構(gòu)，先對該模型結(jié)構(gòu)進行相關(guān)性檢驗并簡化模型結(jié)構(gòu)，直至結(jié)構(gòu)矩陣的所有列之間互不相關(guān)；然后，對簡化模型進行顯著性檢驗并在此基礎(chǔ)上進一步簡化模型結(jié)構(gòu)，直至最終的模型結(jié)構(gòu)為列滿秩，且對應(yīng)的系數(shù)都顯著。

具體步驟為：

(1) 求出結(jié)構(gòu)矩陣的基向量，并以所有的基向量構(gòu)成基矩陣；

(2) 對結(jié)構(gòu)矩陣中除了基矩陣之外的列向量，求出用基矩陣表示的系數(shù)，這些系數(shù)用于構(gòu)建基矩陣對應(yīng)的新的待估計狀態(tài)變量；在求解過程中采用顯著性檢驗，把不顯著的系數(shù)直接置為零；

(3) 根據(jù)基矩陣和觀測向量，求解出新的狀態(tài)變量的值；在求解過程中采用顯著性檢驗，把不顯著的狀態(tài)變量直接置為零。

5 采用高階誤差模型的工具誤差系數(shù)分離

遙外測速度誤差如圖2，可以看出，速度誤差隨著時間積累。

圖2 遙外測速度誤差的擬合結(jié)果Fig.2 Fitting result of remote measurement velocity error

先對高階誤差模型進行相關(guān)性檢驗[9-10]，在環(huán)境函數(shù)矩陣中選取46 個列作為基，這46 個列對應(yīng)的誤差系數(shù)分別為ΔT、k0x、δKax、kyx、Kxxy、Kxyz、K3x、k0y、δKay、kxy、kzy、δK2y、Kyxy、Kyyz、K3y、k0z、δkz、δKaz、kyz、Kzxy、Kzyz、DFx、DFy、DFz、D1x、D1y、D1z、D2x、D2y、D2z、D3y、D3z、D4x、D4z、D5x、D5y、D5z、D6x、D6z、D7x、D7y、D7z、D8x、D8y、D8z、D9z。由這46 個基構(gòu)成基矩陣C′，另外18 個系數(shù)δkx、kzx、K2x、δK2x、Kxxz、δky、K2y、Kyxz、kxz、K2z、δK2z、Kzxz、K3z、D3x、D4y、D6y、D9x、D9y對應(yīng)的列向量都可表示為基的線性組合，這些列向量構(gòu)成矩陣C″。

C″中的每一列都可用C′中基的線性組合表示，設(shè)m=64、d=46，用矩陣方程描述為

式中，p1=3、p2=6、p3=7、p4=8、p5=10、p6=14、p7=18、p8=21、p9=27、p10=29、p11=30、p12=32、p13=34、p14=44、p15=48、p16=54、p17=62、p18=63；q1=1、q2=2、q3=4、q4=5、q5=9、q6=11、q7=12、q8=13、q9=15、q10=16、q11=17、q12=19、q13=20、q14=22、q15=23、q16=24、q17=25、q18=26、q19=28、q20=31、q21=33、q22=35、q23=36、q24=37、q25=38、q26=39、q27=40、q28=41、q29=42、q30=43、q31=45、q32=46、q33=47、q34=49、q35=50、q36=51、q37=52、q38=53、q39=55、q40=56、q41=57、q42=58、q43=59、q44=60、q45=61、q46=64。

矩陣R中各元素的值可通過采用最小二乘法求得，但需強調(diào)的是，經(jīng)過顯著性檢驗后有些元素值可置為零，從而簡化方程。比如，

圖3 δkx 對應(yīng)的列向量元素及其擬合殘差Fig.3 Contrast of δkxvalue and its’ fitting residual error in different conditions

圖4 δky 對應(yīng)的列向量元素及其擬合殘差Fig.4 Contrast of δkyvalue and its’ fitting residual error in different conditions

在求解出矩陣R中的各元素后，考慮到相關(guān)性后的基于遙外測速度差的因果關(guān)系如圖5 所示。

圖5 結(jié)構(gòu)矩陣列相關(guān)時的遙外測速度差因果關(guān)系圖Fig.5 The causality diagram of remote measurement velocity error when the structure matrix columns are correlated

重新定義一組新的誤差系數(shù)：

基于新的誤差系數(shù)的遙外測速度差因果關(guān)系如圖6 所示，其結(jié)構(gòu)矩陣為列滿秩，用線性函數(shù)進行描述見式(28)。

圖6 結(jié)構(gòu)矩陣列滿秩時的遙外測速度差因果關(guān)系圖Fig.6 The causality diagram of remote measurement velocity error based on new error coefficient

再采用最小二乘法并經(jīng)顯著性檢驗后，保留顯著項后的遙外測速度差因果關(guān)系如圖7 所示，表達式為：

圖7 保留顯著項后的遙外測速度差因果關(guān)系圖Fig.7 The causality diagram of remote velocity error retaining significance

最終得到顯著的慣性器件誤差模型為：

除上述系數(shù)顯著外，其余各項誤差系數(shù)都不顯著。誤差系數(shù)不顯著的原因分為三種，第一種原因是在地面上零次項和一次項誤差系數(shù)都已進行了補償，第二種原因是有些高階誤差項在飛行條件下不能充分激勵，第三種原因部分誤差只是數(shù)學(xué)描述而并不存在明確的物理意義。

圖8 遙外測速度誤差補償殘差曲線Fig.8 Residual curve of remote velocity error compensation

6 結(jié)論

本文給出了一種結(jié)合相關(guān)性檢驗和顯著性檢驗的提高慣性制導(dǎo)精度的方法，通過引進環(huán)境函數(shù)矩陣到相關(guān)檢驗中，把相互之間關(guān)聯(lián)的制導(dǎo)工具誤差系數(shù)進行整合，重構(gòu)出新的遙外測速度誤差模型。整合后的新系數(shù)對應(yīng)的環(huán)境函數(shù)矩陣為列滿秩，從而可精確求解誤差系數(shù)，克服了主成分估計、嶺估計等方法不能精確求解的缺點。把顯著性檢驗引入制導(dǎo)工具誤差系數(shù)分離過程中，有利于簡化模型，也有利于分析各制導(dǎo)工具誤差系數(shù)與測量值之間的本質(zhì)特性，并且能夠大幅降低模型的維數(shù)，克服了結(jié)構(gòu)矩陣維數(shù)過多的問題，具有簡單快捷、容易實現(xiàn)的優(yōu)點。將采用本文提出的方法分離得到的誤差系數(shù)代入遙外測速度誤差模型，計算得到的擬合殘差大部分在±0.2 m/s 之內(nèi)，計算結(jié)果表明本文所提出的方法可以有效、準(zhǔn)確地分離慣性制導(dǎo)工具誤差系數(shù)，在慣性制導(dǎo)精度提升領(lǐng)域具有很好的應(yīng)用價值。