淡 鵬，王 丹

(1.宇航動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710043；2.西安衛(wèi)星測(cè)控中心，陜西 西安 710043)

基于多項(xiàng)式插值的外測(cè)數(shù)據(jù)拼幀及平滑問題研究

淡 鵬1，2，王 丹1

(1.宇航動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710043；2.西安衛(wèi)星測(cè)控中心，陜西 西安 710043)

針對(duì)微波統(tǒng)一測(cè)控系統(tǒng)外測(cè)測(cè)量數(shù)據(jù)的插值對(duì)齊、拼幀與平滑計(jì)算問題，結(jié)合工程實(shí)際情況，探討了基于多項(xiàng)式模型的計(jì)算方法，給出了多項(xiàng)式的選取建議，以及測(cè)角等類型測(cè)元數(shù)據(jù)處理的較高精度計(jì)算方法。仿真計(jì)算結(jié)果表明提出的方法可行，不同的插值方法、多項(xiàng)式選取和角度處理方法等對(duì)精度有較大影響，在具體的計(jì)算中需要綜合考慮精度及實(shí)時(shí)性要求的限制。研究結(jié)果對(duì)外測(cè)測(cè)量數(shù)據(jù)的處理有一定的參考價(jià)值。

統(tǒng)一測(cè)控系統(tǒng)；插值；測(cè)距；測(cè)角

0 引言

航天器跟蹤測(cè)量中，微波[1]統(tǒng)一測(cè)控系統(tǒng)[2]外測(cè)[3]測(cè)量數(shù)據(jù)是一類重要的觀測(cè)數(shù)據(jù)，在衛(wèi)星軌道、火箭彈道和返回艙彈道等的計(jì)算中有著重要價(jià)值。但由于測(cè)量機(jī)理及設(shè)備采集方法等方面的原因，此類設(shè)備的外測(cè)測(cè)距[4]、測(cè)速[5]和測(cè)角[6](包括方位角和仰角)等測(cè)元信息在實(shí)際采集后的傳輸中常常不在一幀數(shù)據(jù)中存在，從而給需要同時(shí)使用其中幾類數(shù)據(jù)的計(jì)算帶來困難，如某些實(shí)時(shí)彈道[7]計(jì)算方法需要使用測(cè)距、測(cè)角均有的數(shù)據(jù)，這時(shí)就需要對(duì)其中幾類或全部測(cè)元數(shù)據(jù)進(jìn)行拼幀操作，使之成為測(cè)元“完備”[8]的觀測(cè)數(shù)據(jù)。但這些測(cè)元的時(shí)標(biāo)[9]在修正后或修正前就可能不一致，從而需要對(duì)其進(jìn)行基于時(shí)標(biāo)的插值操作。另外，在數(shù)據(jù)采集頻率較大時(shí)(如20點(diǎn)/s)，為減少數(shù)據(jù)量，有時(shí)還需對(duì)多點(diǎn)數(shù)據(jù)進(jìn)行平滑[10]處理，這些問題從本質(zhì)上來說涉及的是對(duì)這些測(cè)元數(shù)據(jù)的插值及擬合問題。

目前對(duì)外測(cè)數(shù)據(jù)不同測(cè)元進(jìn)行拼幀及平滑問題進(jìn)行深入研究的文獻(xiàn)較難找到，大多是直接使用其數(shù)據(jù)的應(yīng)用算法(如文獻(xiàn)[4-5，11-12])。為此，本文從工程實(shí)際情況出發(fā)，在充分考慮工程計(jì)算面臨的實(shí)際困難，以及精度和可靠性要求的情況下，采用多項(xiàng)式[13]插值的方法對(duì)外測(cè)數(shù)據(jù)拼幀與平滑問題進(jìn)行了分析。

1 外測(cè)數(shù)據(jù)多項(xiàng)式插值及平滑算法

對(duì)外測(cè)數(shù)據(jù)的拼幀及平滑問題，通過實(shí)際的工程計(jì)算可以發(fā)現(xiàn)，不合理的插值算法可能使得拼接后的數(shù)據(jù)計(jì)算的彈道偏差增大。特別是在測(cè)站對(duì)航天器的跟蹤仰角出現(xiàn)最大值的前后(稱為“過頂”，如圖1所示)，此時(shí)角度測(cè)量量變化快、精度低的插值計(jì)算會(huì)使得角度誤差增大，進(jìn)而使得單點(diǎn)測(cè)距與角度數(shù)據(jù)轉(zhuǎn)換為位置時(shí)誤差變大。同時(shí)考慮到在實(shí)時(shí)的外測(cè)數(shù)據(jù)處理中，過于復(fù)雜的插值算法又較難使用。為此具體的實(shí)現(xiàn)中要兼顧可靠性、精度和實(shí)時(shí)性等多個(gè)要求。

圖1 航天器仰角“過頂”示意

1.1 多項(xiàng)式插值及拼幀方法

對(duì)外測(cè)數(shù)據(jù)中不同測(cè)量元素?cái)?shù)據(jù)的插值與拼幀、平滑操作的基本算法是曲線擬合與插值，這方面的方法很多，如較常用的拉格朗日插值[13]、最小二乘曲線擬合插值[14]、樣條插值[13]和拋物線插值[13]等。其中，拉格朗日插值具有算法簡單、實(shí)現(xiàn)容易的特點(diǎn)，而最小二乘[15]曲線擬合插值在消除野值方面具有較好表現(xiàn)。為此本文在多項(xiàng)式插值與擬合計(jì)算中，主要采用了拉格朗日與最小二乘曲線擬合2種方法，相關(guān)方法實(shí)現(xiàn)可參考文獻(xiàn)[13]。

當(dāng)測(cè)距、測(cè)角和測(cè)速等測(cè)元數(shù)據(jù)插值到同一時(shí)間點(diǎn)后，即可按時(shí)間基準(zhǔn)進(jìn)行拼幀操作。在插值及拼幀操作中需要注意以下方面：

① 時(shí)標(biāo)的跨天處理。在外測(cè)測(cè)元數(shù)列插值中，需要解決的一個(gè)問題是時(shí)標(biāo)的跨天問題。處理方法為：對(duì)一時(shí)間點(diǎn)為T1,T2,…,Tn的測(cè)量序列，插值計(jì)算TK處的值，可將其時(shí)間點(diǎn)轉(zhuǎn)換為T1-TK,T2-TK,…,Tn-TK，這樣只需插值T=0處結(jié)果。此處通過對(duì)時(shí)間做差，避免了跨天問題。

② 插值中的時(shí)標(biāo)對(duì)齊策略。在不同測(cè)元數(shù)據(jù)對(duì)齊時(shí)，可采用2種插值策略：一是對(duì)每一個(gè)測(cè)量量都進(jìn)行插值，然后以時(shí)間驅(qū)動(dòng)查找相同的時(shí)標(biāo)點(diǎn)進(jìn)行拼接；二是以某一個(gè)測(cè)量量時(shí)標(biāo)為基準(zhǔn)，將其他測(cè)量量插值到此時(shí)標(biāo)處，這樣可減小某一維數(shù)據(jù)的插值誤差。如可采用以測(cè)距數(shù)據(jù)時(shí)標(biāo)為基準(zhǔn)，將測(cè)角、測(cè)速插值到其時(shí)標(biāo)處然后進(jìn)行拼接。

1.2 計(jì)算中的角度數(shù)據(jù)處理

角度數(shù)據(jù)是外測(cè)數(shù)據(jù)插值與平滑計(jì)算中需要特別對(duì)待的一類。此類數(shù)據(jù)在數(shù)值上可能是不連續(xù)的，如方位角的“過北”現(xiàn)象(從360°變?yōu)?°或相反過程)，且插值精度低時(shí)會(huì)使誤差變大。為此給出下面的轉(zhuǎn)換插值算法。

① 對(duì)仰角E：若直接對(duì)角度進(jìn)行插值擬合則偏差太大，故改為E的函數(shù)f(E)進(jìn)行插值計(jì)算(f(E)通常取為正弦函數(shù))，最后再解算出E。

此時(shí)，函數(shù)f(E)變化率相對(duì)較小，一定程度上可減小插值誤差。

② 對(duì)方位角A：在方位角“過北”情況下(即方位角從第1像限進(jìn)入第4像限、或者從第4像限進(jìn)入第1像限時(shí))，插值計(jì)算時(shí)需要進(jìn)行角度連續(xù)性判斷，判斷方法為，如果前后兩點(diǎn)方位角滿足|Ak-1-Ak|>180.0，則認(rèn)為過北。

為簡化計(jì)算并正確處理過北情況，對(duì)A的處理可采用以下2種方法：

方法1：采用三角函數(shù)法進(jìn)行插值計(jì)算，具體為對(duì)正弦與余弦函數(shù)sin(A)、cos(A)分別進(jìn)行插值計(jì)算，然后再借用反三角函數(shù)arctan2解算出A。

關(guān)于角度與三角函數(shù)插值的誤差問題，通過對(duì)某任務(wù)某觀測(cè)站過頂前后的實(shí)測(cè)跟蹤數(shù)據(jù)分別采用角度直接插值以及三角函數(shù)插值2種方法計(jì)算，然后做差，得到的差值曲線如圖2所示。

圖2 測(cè)角數(shù)據(jù)2種方法插值結(jié)果互差散點(diǎn)

從圖2中可看出，此樣本數(shù)據(jù)2種方法計(jì)算結(jié)果，方位角最大差值為0.000 115°，仰角最大差值為0.000 283°，并且仰角偏差大的地方出現(xiàn)在其變化率大的地方，2種方法計(jì)算結(jié)果差別很小。

考慮到航天器飛行過程是位置、速度連續(xù)變化的過程，對(duì)TK處外測(cè)觀測(cè)量的插值最終要?dú)w結(jié)為TK處飛行狀態(tài)量的插值。因此也可從二者之間的轉(zhuǎn)換關(guān)系比較直接使用角度與使用三角函數(shù)時(shí)的精度。

對(duì)測(cè)站地平坐標(biāo)系[4]的球坐標(biāo)觀測(cè)量ρk,Ak,Ek和目標(biāo)在地心固連坐標(biāo)系[5]的坐標(biāo)可表示為

(1)

式中，M為測(cè)站地平系到地心固連系的轉(zhuǎn)換矩陣；R為測(cè)站在地心固連系下的坐標(biāo)。當(dāng)測(cè)站確定時(shí)，M與R為確定量，而測(cè)距ρk的插值結(jié)果一定情況下，位置的插值精度直接決定于測(cè)角三角函數(shù)的插值精度。

1.3 多項(xiàng)式平滑算法

2 仿真計(jì)算與分析

為驗(yàn)證算法有效性及對(duì)不同條件下的計(jì)算效果進(jìn)行對(duì)比，給定一組軌道初值，外推計(jì)算出某測(cè)站在過頂前后非整秒點(diǎn)的理論測(cè)量值進(jìn)行曲線擬合插值計(jì)算。

考慮到在實(shí)時(shí)工程計(jì)算中，對(duì)數(shù)據(jù)處理時(shí)延要求嚴(yán)格，為此在時(shí)標(biāo)插值計(jì)算中最多使用插值點(diǎn)后1 s內(nèi)的數(shù)據(jù)進(jìn)行計(jì)算。

2.1 拉格朗日插值測(cè)試

仿真數(shù)據(jù)為1 點(diǎn)/s，作以下2種計(jì)算：

算例1：采用拉格朗日插值方法，使用了插值點(diǎn)后的1點(diǎn)數(shù)據(jù)(即進(jìn)行內(nèi)插)，分別計(jì)算了總點(diǎn)數(shù)為2點(diǎn)、3點(diǎn)和4點(diǎn)擬合結(jié)果，并與理論測(cè)量值做差。

算例2：只使用插值點(diǎn)時(shí)間之前的數(shù)據(jù)(即進(jìn)行外插)，外插后1 s的數(shù)據(jù)，分別計(jì)算了總點(diǎn)數(shù)為2點(diǎn)、3點(diǎn)和4點(diǎn)計(jì)算結(jié)果，并與理論測(cè)量值做差。

計(jì)算結(jié)果如圖3所示。需要說明的是，外插時(shí)2點(diǎn)計(jì)算結(jié)果中某些值超出范圍，未能在圖中顯示出來。

從圖3中可以看出，對(duì)于內(nèi)插，3點(diǎn)與4點(diǎn)插值誤差明顯小于2點(diǎn)插值結(jié)果，而4點(diǎn)插值結(jié)果誤差已經(jīng)很小，若采用更多點(diǎn)插值，誤差并沒有太多提高，相反使用過多點(diǎn)還可能使結(jié)果出現(xiàn)病態(tài)。而對(duì)于外插，4點(diǎn)插值結(jié)果要好于3點(diǎn)，而3點(diǎn)又明顯好于2點(diǎn)結(jié)果?？傮w來說，內(nèi)插好于外插。

圖3 拉格朗日插值對(duì)比

2.2 最小二乘擬合插值測(cè)試

仿真1 點(diǎn)/s的數(shù)據(jù)，分別使用10點(diǎn)、7點(diǎn)數(shù)據(jù)做2次和3次多項(xiàng)式最小二乘擬合插值(內(nèi)插前1 s的數(shù)據(jù))，然后將插值結(jié)果與理論值做差，所得差值曲線如圖4所示。

圖4 最小二乘插值對(duì)比

從圖4中可以看出，3次多項(xiàng)式最小二乘插值效果要好于2次多項(xiàng)式。而10點(diǎn)數(shù)據(jù)最小二乘插值結(jié)果并不比7點(diǎn)好，這是因?yàn)辄c(diǎn)數(shù)太多時(shí)，使用的插值點(diǎn)時(shí)標(biāo)之前的數(shù)據(jù)過多(出于實(shí)時(shí)計(jì)算時(shí)延考慮，此處只使用了目標(biāo)插值點(diǎn)后1 s的數(shù)據(jù)，其余均為之前的數(shù)據(jù))，使得擬合值對(duì)前面點(diǎn)的方差更小，曲線更趨近于前面的點(diǎn)。

2.3 多項(xiàng)式最小二乘平滑測(cè)試

仿真20 點(diǎn)/s數(shù)據(jù)，采用每2 s數(shù)據(jù)進(jìn)行多項(xiàng)式中心最小二乘平滑計(jì)算，平滑后數(shù)據(jù)與理論計(jì)算值做差，可得到測(cè)距差在1 cm范圍內(nèi)，測(cè)角差在1″范圍內(nèi)，滿足計(jì)算的精度要求。

2.4 結(jié)論

從上面的仿真計(jì)算可得出如下結(jié)論：

① 本文給出的基于多項(xiàng)式的外測(cè)數(shù)據(jù)插值與平滑算法是有效的，通過三角函數(shù)來實(shí)現(xiàn)角度數(shù)據(jù)的插值，可有效避免直接使用角度造成的不連續(xù)問題；

② 在外測(cè)數(shù)據(jù)插值計(jì)算時(shí)，推薦使用3次多項(xiàng)式插值，當(dāng)采用拉格朗日插值算法時(shí)，可使用4點(diǎn)公式計(jì)算(即3次多項(xiàng)式)；

③ 多項(xiàng)式最小二乘法的插值對(duì)齊與平滑算法對(duì)野值的剔除有一定效果，但使用時(shí)應(yīng)盡量進(jìn)行內(nèi)插，外插時(shí)精度降低較大。

同時(shí)應(yīng)該看到，擬合結(jié)果除了與使用點(diǎn)數(shù)、擬合方法有關(guān)外，還與數(shù)據(jù)的變化率和插值位置等因素有關(guān)，因此上面給出的結(jié)果只是針對(duì)于給定的計(jì)算條件。

綜合對(duì)比圖3和圖4可以發(fā)現(xiàn)，最小二乘曲線擬合插值結(jié)果誤差并不一定比拉格朗日結(jié)果好，在實(shí)時(shí)外測(cè)數(shù)據(jù)處理中，使用拉格朗日插值已經(jīng)能夠滿足精度要求。但是考慮到觀測(cè)值的野值問題，在數(shù)據(jù)質(zhì)量不好的情況下使用最小二乘擬合方法結(jié)果可能會(huì)好一些。

3 結(jié)束語

本文對(duì)基于多項(xiàng)式算法的外測(cè)數(shù)據(jù)插值、拼幀及平滑方法進(jìn)行了探討，對(duì)計(jì)算中的有關(guān)細(xì)節(jié)問題進(jìn)行了分析，計(jì)算結(jié)果表明了本文算法的有效性。

應(yīng)該看到，在實(shí)際工程計(jì)算中，外測(cè)數(shù)據(jù)實(shí)時(shí)拼幀處理中的插值對(duì)齊問題需要綜合考慮時(shí)延要求和插值位置等因素。綜合上面的計(jì)算結(jié)果，認(rèn)為實(shí)時(shí)插值時(shí)使用的點(diǎn)數(shù)不宜太多，應(yīng)盡量避免外插的發(fā)生(但會(huì)增大計(jì)算延遲)，同時(shí)簡單的拉格朗日插值在精度上可以滿足與最小二乘相當(dāng)?shù)木?數(shù)據(jù)質(zhì)量不好時(shí)差些，但這可通過在具體的數(shù)據(jù)使用過程中，應(yīng)用濾波等方法進(jìn)行二次處理來解決)。下一步將重點(diǎn)在測(cè)量數(shù)據(jù)的濾波處理等方面進(jìn)行研究。

[1] ZAKA I,REHMAN H U,SHAH S I.PSO- and GA-based Jammer Excision in CDMA[J].Journal of Circuits Systems and Computers,2010,19(1):123-138.

[2] 趙楠.統(tǒng)一擴(kuò)頻測(cè)控系統(tǒng)中關(guān)鍵技術(shù)研究[D].哈爾濱：哈爾濱工業(yè)大學(xué)，2011.

[3] 茅永興,馬靜遠(yuǎn),鐘德安,等.基于速度增量修正的衛(wèi)星入軌段外測(cè)數(shù)據(jù)定軌新方法[J].宇航學(xué)報(bào)，2016，37(9)：1 049-1 055.

[4] 淡鵬,李恒年,李志軍.應(yīng)用三向測(cè)量數(shù)據(jù)的深空探測(cè)器實(shí)時(shí)濾波定位算法[J].航天器工程,2015,24(2):21-26.

[5] 淡鵬,李恒年,張智斌.一種火箭外測(cè)彈道實(shí)時(shí)重建的自適應(yīng)濾波算法[J].彈箭與制導(dǎo)學(xué)報(bào),2013,33(6)：188-190.

[6] 王英玲.低仰角跟蹤技術(shù)研究及應(yīng)用[J].無線電工程,2008,38(9):43-46.

[7] 胡亞男,甘友誼,王盛璽.衛(wèi)星發(fā)射任務(wù)主動(dòng)段實(shí)時(shí)彈道算法研究[J].無線電工程,2013,43(4):28-31.

[8] 淡鵬,李恒年,張定波,等.基于多元非完備信息的實(shí)時(shí)濾波定軌方法[J].飛行力學(xué),2014,32(3)：283-288.

[9] 王金華,劉魁星,韋欣榮.基于矩陣構(gòu)造的偽碼測(cè)距時(shí)標(biāo)偏差標(biāo)定方法[J].無線電工程,2016,46(12):47-50.[10] 梁健,王京京.基于改進(jìn)Kalman濾波算法的航跡平滑[J].無線電工程,2015,45(5):20-23.

[11] 張艷,郭軍海,慈穎.基于動(dòng)力學(xué)約束的實(shí)時(shí)彈道滑動(dòng)處理方法[J].飛行器測(cè)控學(xué)報(bào),2011,30(3):56-60.

[12] 雷鳴,張玉如.自適應(yīng)濾波方法在外彈道實(shí)時(shí)測(cè)量中的應(yīng)用[J].艦船電子工程,2002(3):18-24.

[13] 石瑞民.數(shù)值計(jì)算[M].北京：高等教育出版社，2004.

[14] 王易南,閆杰,溫琦.基于最小二乘的航天器天文導(dǎo)航定姿算法研究[J].電子設(shè)計(jì)工程,2014,23(5):58-59.

[15] 米芳彬,徐小剛,梁健.基于有限元記憶最小二乘的雷達(dá)誤差配準(zhǔn)算法[J].無線電工程,2016,46(2):65-68.

[16] 劉玉明,張寶華.一類FIR濾波器的加權(quán)最小方差設(shè)計(jì)[J].無線電工程,2000,30(4):52-54.

Discussion on Tracking Data Completing & Smoothing Using Polynomial Interpolation

DAN Peng1,2，WANG Dan1

(1.StateKeyLaboratoryofAstronauticDynamics,Xi’anShaanxi710043,China;2.Xi’anSatelliteControlCenter,Xi’anShaanxi710043,China)

For tracking data interpolation,completing and smoothing of microwave unified TT&C system,a polynomial calculation method is discussed,some suggestions for the selection of polynomial formula as well as some high-precision processing method for bearing measurement are provided.The simulation results show that the method given in this paper is feasible,and different formulas,interpolation methods and bearing processing methods will result in apparent accuracy deviation.In practical engineering calculation,the result precision and real-time requirement,etc.need to be considered in the processing of tracking data.The results provide a reference for the processing of external tracking data.

unified TT&C system;interpolation;ranging;bearing measurement

10.3969/j.issn.1003-3106.2017.07.05

淡鵬,王丹.基于多項(xiàng)式插值的外測(cè)數(shù)據(jù)拼幀及平滑問題研究[J].無線電工程,2017,47(7):20-24.[DAN Peng,WANG Dan.Discussion on Tracking Data Completing & Smoothing Using Polynomial Interpolation[J].Radio Engineering,2017,47(7):20-24.]

2017-01-10

國家自然科學(xué)基金資助項(xiàng)目(61302098)。

V412.4

A

1003-3106(2017)07-0020-05

淡 鵬 男，(1979—)，碩士，工程師。主要研究方向：航天器數(shù)據(jù)處理技術(shù)。

王 丹 女，(1980—)，碩士，工程師。主要研究方向：航天器數(shù)據(jù)處理技術(shù)。