陳晶，李國(guó)友，黃維，張新蘭

(中國(guó)人民解放軍95948部隊(duì)，甘肅 蘭州　732750)

?

含特征點(diǎn)飛行器位置參數(shù)的樣條微分平滑方法

陳晶，李國(guó)友，黃維，張新蘭

(中國(guó)人民解放軍95948部隊(duì)，甘肅 蘭州732750)

為了獲取準(zhǔn)確的飛行器速度參數(shù)，根據(jù)含特征點(diǎn)的飛行器軌道特點(diǎn)，提出了一種基于最優(yōu)節(jié)點(diǎn)樣條逼近的微分平滑算法。首先給出了含特征點(diǎn)軌道位置參數(shù)樣條微分平滑問(wèn)題的數(shù)學(xué)模型；其次介紹了特征點(diǎn)的識(shí)別方法；然后建立了基于混沌粒子群算法的樣條微分平滑算法；最后給出了典型情況下的仿真結(jié)果。仿真結(jié)果表明：該方法可以有效減小隨機(jī)誤差和截?cái)嗾`差的影響，提高含特征點(diǎn)飛行器軌道速度參數(shù)計(jì)算精度。

飛行器；位置參數(shù)；微分平滑；特征點(diǎn)；樣條節(jié)點(diǎn)；粒子群優(yōu)化

0　引言

在光測(cè)數(shù)據(jù)事后處理中，通常采用多臺(tái)光學(xué)經(jīng)緯儀進(jìn)行空間交會(huì)獲取飛行器位置參數(shù)，再進(jìn)行正交多項(xiàng)式中心微分平滑，獲取飛行器速度參數(shù)。選擇多項(xiàng)式階數(shù)和平滑區(qū)間長(zhǎng)度需要綜合考慮隨機(jī)誤差和截?cái)嗾`差影響，在運(yùn)動(dòng)變化較大時(shí)，特別是在含有特征點(diǎn)(如主動(dòng)段關(guān)機(jī)點(diǎn)、級(jí)間分離點(diǎn))的段落，采用多項(xiàng)式微分平滑存在較大的截?cái)嗾`差，影響速度參數(shù)的準(zhǔn)確性[1]。

文獻(xiàn)[2-7]的研究表明，飛行器的軌道位置參數(shù)可用3次B樣條函數(shù)高精度逼近。利用樣條函數(shù)的可微性可以獲取飛行器速度參數(shù)。文獻(xiàn)[8-10]通過(guò)優(yōu)化樣條節(jié)點(diǎn)分布適應(yīng)目標(biāo)運(yùn)動(dòng)規(guī)律，在運(yùn)動(dòng)變化較大段落節(jié)點(diǎn)分布較密，減小了截?cái)嗾`差；在運(yùn)動(dòng)變化緩慢段落節(jié)點(diǎn)分布稀疏，減小了隨機(jī)誤差。但在樣條節(jié)點(diǎn)優(yōu)化過(guò)程中，需要利用對(duì)數(shù)變換將有約束問(wèn)題轉(zhuǎn)化為無(wú)約束問(wèn)題，并需要計(jì)算待優(yōu)化函數(shù)的雅克比矩陣，算法復(fù)雜，并且沒(méi)有考慮飛行器軌道特征點(diǎn)識(shí)別。

本文針對(duì)含有特征點(diǎn)的飛行器軌道位置參數(shù)的樣條微分平滑問(wèn)題，采用遙測(cè)數(shù)據(jù)或光測(cè)圖像確定飛行器軌道特征點(diǎn)時(shí)刻，將特征點(diǎn)作為固定樣條節(jié)點(diǎn)，采用混沌粒子群算法對(duì)其余節(jié)點(diǎn)進(jìn)行優(yōu)化分布，實(shí)現(xiàn)了對(duì)含特征點(diǎn)飛行器軌道位置參數(shù)的高精度樣條逼近，減小了截?cái)嗾`差和隨機(jī)誤差影響，并且不需要計(jì)算優(yōu)化函數(shù)的雅克比矩陣，降低了算法復(fù)雜程度。

1　問(wèn)題的數(shù)學(xué)描述

飛行器軌道受力變化導(dǎo)致加速度突變的點(diǎn)被稱(chēng)為軌道特征點(diǎn)，如主動(dòng)段關(guān)機(jī)點(diǎn)、級(jí)間分離點(diǎn)等，當(dāng)特征點(diǎn)被確定為樣條節(jié)點(diǎn)時(shí)，飛行器軌道空間位置參數(shù)可以用3次B樣條函數(shù)進(jìn)行高精度表示[11]。

假設(shè)多臺(tái)光電經(jīng)緯儀在t1～tm時(shí)間段內(nèi)對(duì)飛行器進(jìn)行跟蹤測(cè)量，進(jìn)行外測(cè)數(shù)據(jù)事后處理通過(guò)空間交會(huì)得到軌道位置參數(shù)(x(ti),y(ti),z(ti))T,ti∈[t1,tm]。

下面以y方向位置參數(shù)的樣條微分平滑為例進(jìn)行說(shuō)明。飛行器ti時(shí)刻y方向位置參數(shù)y(ti)可以表示為

(1)

式中：Bj(τ,ti)為以τ為樣條內(nèi)節(jié)點(diǎn)的3次規(guī)范B樣條基函數(shù)

(2)

βj為軌道參數(shù)y方向的樣條系數(shù)；η(ti)為飛行器軌道y方向位置參數(shù)的測(cè)量誤差。

定義A(τ)為以τ為樣條內(nèi)節(jié)點(diǎn)的樣條基矩陣

(3)

則軌道y方向位置參數(shù)向量可以表示為

Y=A(τ)β+η，

(4)

式中：Y=(Y(t1),…,Y(tm))T為外測(cè)數(shù)據(jù)事后處理得到的飛行器y方向位置參數(shù)向量；β=(β1,…，β3M)T軌道y方向位置參數(shù)樣條系數(shù)向量；η=(η(t1),…,η(tm))T為y方向位置參數(shù)的誤差向量。

通過(guò)交會(huì)計(jì)算得到的飛行器位置參數(shù)的誤差并不是等方差的，定義

E(ηηT)=P，P=diag(σ2(t1),…,σ2(tm)),

式中：σ2(ti)為ti時(shí)刻飛行器y方向位置參數(shù)隨機(jī)誤差方差。

定義待優(yōu)化目標(biāo)函數(shù)為殘差平方和

Q(τ,β)=‖P-1/2Y-P-1/2A(τ)β‖2，

(5)

位置參數(shù)的樣條微分平滑問(wèn)題就是對(duì)τ，β進(jìn)行估計(jì)，使得

(6)

(7)

利用樣條函數(shù)的可微性可以對(duì)軌道位置參數(shù)進(jìn)行微分計(jì)算得到速度參數(shù)

(8)

含特征點(diǎn)飛行器位置參數(shù)的樣條微分平滑問(wèn)題等價(jià)于先對(duì)樣條內(nèi)節(jié)點(diǎn)序列τ進(jìn)行優(yōu)化，然后對(duì)軌道參數(shù)的樣條系數(shù)向量β進(jìn)行估計(jì)的優(yōu)化問(wèn)題。在優(yōu)化過(guò)程中，不僅需要考慮樣條內(nèi)節(jié)點(diǎn)序列τ自身的時(shí)間先后順序約束條件，即t1

2　軌道參數(shù)特征點(diǎn)的確定方法

實(shí)際飛行器外測(cè)數(shù)據(jù)事后處理過(guò)程中，軌道特征點(diǎn)序列τ*并非已知量，準(zhǔn)確定位飛行器軌道特征點(diǎn)成為軌道位置參數(shù)樣條微分平滑的前提。

飛行器的主動(dòng)段關(guān)機(jī)、級(jí)間分離、二級(jí)點(diǎn)火等特征時(shí)刻，飛行器光測(cè)圖像會(huì)有明顯變化，可以直接通過(guò)光測(cè)圖像進(jìn)行特征點(diǎn)確定。當(dāng)飛行器具有遙測(cè)數(shù)據(jù)時(shí)，可以直接利用遙測(cè)的視加速度數(shù)據(jù)確定軌道特征點(diǎn)[12]。通過(guò)飛行器平臺(tái)安裝的3個(gè)加速度表，測(cè)量沿慣性坐標(biāo)系3個(gè)方向的視加速度，視加速度發(fā)生突變的點(diǎn)就是飛行器軌道的特征點(diǎn)。值得注意的是，由于遙測(cè)系統(tǒng)和外測(cè)系統(tǒng)在采樣頻率、時(shí)間零點(diǎn)等方面存在差異，需要對(duì)遙測(cè)數(shù)據(jù)進(jìn)行必要的對(duì)齊和修正。

3　樣條微分平滑的粒子群優(yōu)化算法

特征點(diǎn)序列確定后的位置參數(shù)樣條微分平滑問(wèn)題是一個(gè)帶有復(fù)雜約束條件的非線性?xún)?yōu)化問(wèn)題，傳統(tǒng)的高斯-牛頓法難以完成尋優(yōu)計(jì)算[13]?；煦缌Ｗ尤核惴ㄊ且环N群體智能優(yōu)化算法，通過(guò)引入混沌思想，既具有粒子群算法簡(jiǎn)單、易于編程實(shí)現(xiàn)的優(yōu)點(diǎn)，又克服了復(fù)雜函數(shù)優(yōu)化的早熟收斂問(wèn)題，可用于解決特征點(diǎn)序列確定后樣條節(jié)點(diǎn)優(yōu)化這類(lèi)復(fù)雜的非線性?xún)?yōu)化問(wèn)題[14]。

3.1適應(yīng)度計(jì)算

在粒子群算法中，需要利用適應(yīng)度值判斷粒子的優(yōu)劣，對(duì)于特征點(diǎn)序列確定后的樣條微分平滑問(wèn)題，采用樣條逼近殘差平方和Q(τ′,β)作為適應(yīng)度值。

計(jì)算粒子Γi的適應(yīng)度值步驟如下：

(1) 將Γi和τ*合并，按從小到大排列得到樣條內(nèi)節(jié)點(diǎn)序列為τ=(T5,T6,…,TM)T，進(jìn)行節(jié)點(diǎn)延拓后得到總的節(jié)點(diǎn)序列T1

(2) 計(jì)算樣條基矩陣，根據(jù)式(2)計(jì)算樣條內(nèi)節(jié)點(diǎn)序列τ對(duì)應(yīng)的樣條基函數(shù)，代入式(3)計(jì)算得到樣條基矩陣A(τ)。

(3) 將A(τ)代入式(4)，進(jìn)行馬爾可夫估計(jì)得到樣條系數(shù)向量的估計(jì)值

(9)

(10)

3.2粒子位置更新

粒子i根據(jù)下面的公式更新位置和速度

(11)

(12)

式中：c1和c2為學(xué)習(xí)因子常數(shù)；r1和r2為[0,1]之間的隨機(jī)數(shù)；εd為小于飛行器軌道位置參數(shù)時(shí)間采樣間隔的小量。

3.3全局最優(yōu)位置的混沌優(yōu)化

在復(fù)雜優(yōu)化問(wèn)題中，粒子群算法易產(chǎn)生早熟收斂問(wèn)題，陷入局部極值點(diǎn)，因此引入混沌粒子群算法克服這一問(wèn)題。根據(jù)全局最優(yōu)位置采用混沌映射產(chǎn)生混沌序列種群，然后利用混沌序列種群中的最優(yōu)粒子隨機(jī)替代當(dāng)前種群中的某個(gè)粒子[15]。

(13)

3.4算法流程

采用混沌粒子群算法進(jìn)行特征點(diǎn)確定后的樣條節(jié)點(diǎn)位置優(yōu)化的流程為

步驟1：初始化粒子群算法基本參數(shù)。根據(jù)自由樣條內(nèi)節(jié)點(diǎn)個(gè)數(shù)D，計(jì)算初始樣條節(jié)點(diǎn)間隔h=(tm-t1)/D，設(shè)置最大粒子速度為，vmax=h設(shè)置學(xué)習(xí)因子c1，c2，設(shè)置粒子群種群規(guī)模maxS，設(shè)置最大進(jìn)化代數(shù)maxE，設(shè)置混沌迭代最大次數(shù)maxL。

步驟2：初始化種群，隨機(jī)生成maxS個(gè)自由樣條內(nèi)節(jié)點(diǎn)粒子

Γi=t1+rand(1,D)·(tm-t1)，

(14)

式中：rand(1,D)為取值介于[0,1]之間的1×D維隨機(jī)數(shù)向量。

粒子Γi的速度初始化為

Vi=-vmax+rand(1,D)·2vmax.

(15)

初始化粒子Γi的最好位置為Pi=Γi，計(jì)算每個(gè)粒子對(duì)應(yīng)的適應(yīng)度值Ei，初始化每個(gè)粒子的最優(yōu)適應(yīng)度值Eibest=Ei，初始化全局最優(yōu)適應(yīng)度Ebest為向量(E1best,…,EDbest)T的最小值，對(duì)應(yīng)的粒子存入種群的全局最優(yōu)位置Pg，進(jìn)化代數(shù)設(shè)置為Gn=1。

步驟3：根據(jù)式(11)，(12)對(duì)每個(gè)粒子位置和速度進(jìn)行更新，計(jì)算每個(gè)粒子的適應(yīng)度Ei，比較Ei和Eibest，若Ei

步驟4：對(duì)全局最優(yōu)位置Pg進(jìn)行混沌優(yōu)化，用混沌種群中的最優(yōu)粒子隨機(jī)代替原種群中的某個(gè)粒子，進(jìn)化代數(shù)Gn=Gn+1。

步驟5：判斷是否達(dá)到最大進(jìn)化代數(shù)或全局最優(yōu)適應(yīng)度值是否滿足精度要求，若沒(méi)有滿足則轉(zhuǎn)入步驟3繼續(xù)進(jìn)行迭代，否則轉(zhuǎn)入步驟6。

步驟6：輸出全局最優(yōu)適應(yīng)度Ebest和全局最優(yōu)位置Pg。

4　仿真驗(yàn)證

4.1仿真條件

利用某型飛行器0～60 s的理論軌道數(shù)據(jù)進(jìn)行仿真，其中包含兩處特征點(diǎn)。根據(jù)光電經(jīng)緯儀站址生成2臺(tái)光測(cè)設(shè)備的角度觀測(cè)數(shù)據(jù)，添加小于15″的隨機(jī)誤差。根據(jù)2臺(tái)光電經(jīng)緯儀的仿真角度數(shù)據(jù)進(jìn)行空間交會(huì)定位，得到飛行器軌道位置參數(shù)測(cè)量值(輸出時(shí)間間隔為0.1 s)。

4.2仿真方案

分別采用如下3種方法對(duì)y方向位置參數(shù)測(cè)量值進(jìn)行微分平滑，得到y(tǒng)方向速度分量vy，與y方向速度分量理論值進(jìn)行比較。

方案1：采用21點(diǎn)二階多項(xiàng)式進(jìn)行多項(xiàng)式中心微分平滑；

方案2：采用文獻(xiàn)[2]方法進(jìn)行樣條微分平滑，樣條內(nèi)節(jié)點(diǎn)數(shù)設(shè)置為20個(gè)；

方案3：采用本文方法進(jìn)行樣條微分平滑，樣條內(nèi)節(jié)點(diǎn)數(shù)設(shè)置為20個(gè)(其中包含2個(gè)特征點(diǎn))。

4.3仿真結(jié)果

3種微分平滑方法得到的與理論值之間差別的均方差如表1所示。

表1　vy誤差均方差對(duì)比

特征點(diǎn)時(shí)刻vy與理論值之間差別的絕對(duì)值如表2所示。

表2　特征點(diǎn)時(shí)刻vy誤差絕對(duì)值對(duì)比

為直觀起見(jiàn)，3種方法得到的vy與理論值在特征點(diǎn)1的局部放大效果如圖1所示。

圖1　不同方法在特征點(diǎn)1的速度分量曲線Fig.1　Velocity curve of characteristic point 1 with different methods

由表1，2和圖1仿真結(jié)果可知，與另外2種方法相比，采用本文方法對(duì)飛行器位置參數(shù)進(jìn)行微分平滑得到的速度分量更接近理論值，特別是在特征點(diǎn)時(shí)刻優(yōu)勢(shì)更為明顯，有效減小了截?cái)嗾`差和隨機(jī)誤差。

5　結(jié)束語(yǔ)

本文在含有特征點(diǎn)的飛行器非平穩(wěn)飛行段落，利用遙測(cè)數(shù)據(jù)或者光測(cè)圖像預(yù)先分析出這些特征點(diǎn)，以特征點(diǎn)為樣條內(nèi)節(jié)點(diǎn)，采用粒子群算法對(duì)其余樣條內(nèi)節(jié)點(diǎn)進(jìn)行優(yōu)化分布，可以有效減小軌道位置參數(shù)樣條微分平滑的隨機(jī)誤差和截?cái)嗾`差，可用于光測(cè)數(shù)據(jù)事后處理的飛行器位置參數(shù)微分平滑，提高速度參數(shù)的計(jì)算精度。

[1]柴敏，胡紹林，郭小紅.外測(cè)彈道截?cái)嗾`差修正方法[J].導(dǎo)彈與航天運(yùn)載技術(shù)，2006(5)：53-57.CHAI Min，HU Shao-lin，GUO Xiao-hong. Research on Algorithm for Truncation Error of Exterior Measured Trajectory[J]. Missile and Space Vehcile，2006(5)：53-57.

[2]王菖，段曉君，王正明，等.B樣條分頻理論及在外測(cè)數(shù)據(jù)處理中的應(yīng)用[J].彈道學(xué)報(bào)，2010，22(3)：99-102.

WANG Chang, DUAN Xiao-jun, WANG Zheng-ming, et al. B-Spline Separating Signals Method and Its Application in Measured Data Post-Processing[J]. Jounal of Ballistics, 2010, 22(3):99-102.

[3]周海銀，李董輝，蔣月評(píng). 基于多源信息融合的飛行器高精度定位方法[J]. 湖南大學(xué)學(xué)報(bào)：自然科學(xué)版，2007，34(4)：37-40.

ZHOU Hai-yin，LI Dong-hui，JIANG Yue-ping. High Precision Method for Determining the Position of Aircraft based on Multi-Sensors Information Fusion[J]. Journal of Hunan University：Natural Sciences ed，2007，34(4)：37-40.

[4]張艷，石晟瑋，郭軍海. 基于樣條約束的彈道融合抗差估計(jì)方法[J]. 彈道學(xué)報(bào)，2011，23(1)：40-44.

ZHANG Yan， SHI Sheng-wei， GUO Jun-hai. Robust Estimation for Trajectory Fusion Based on Spline Restraint[J]. Journal of Ballistics，2011，23(1)：40-44.

[5]賀明科，王正明. 光學(xué)接力跟蹤的彈道融合[J]. 彈道學(xué)報(bào)，2002，14(1)：1-7.

HE Ming-ke，WANG Zheng-ming. Fusion of Trajectory Tracked by Relay Optical Measure System[J]. Journal of Ballistics，2002，14(1)：1-7.

[6]朱武宣，高耀文. 基于樣條約束“EMBET”的再入軌道測(cè)量數(shù)據(jù)融合方法[J]. 飛行器測(cè)控學(xué)報(bào)，2005，24(6)：49-53.

ZHU Wu-xuan，GAO Yao-wen. A Data Fusion Method for Reentry Ballistic Measurement with Spline Restraint EMBET[J]. Journal of Spacecraft TT&C Technology，2005，24(6)：49-53.

[7]趙文策，周海銀，段曉君，等. 基于測(cè)角元和高程信息融合的彈道定位求速[J]. 彈道學(xué)報(bào)，2004，16(2)：27-32.

ZHAO Wen-ce，ZHOU Hai-yin，DUAN Xiao-jun，et al. Trajectory Positioning and Velocity Determination Based on Angle Measurement and Height Information Fusion[J]. Journal of Ballistics，2004，16(2)：27-32.

[8]郭文勝，宮志華，董立濤，等.雷達(dá)測(cè)距擬合微分求速方法研究[J].現(xiàn)代雷達(dá)，2010，32(8)：33-37.

GUO Wen-sheng, GONG Zhi-hua, DONG Li-tao, et al. A Study on Getting Radar Radial Velocity by Range Fitting and Differential Method[J]. Modern Radar, 2010, 32(8):33-37.

[9]易東云，朱炬波，王正明.三頻帶信號(hào)分離的樣條函數(shù)方法[J].電子學(xué)報(bào)，1999(1)：54-57.

YI Dong-yun, ZHU Ju-bo, WANG Zheng-ming. Three Bands Seperating Signal Method with Spline Function[J]. Jounal of Electronics, 1999(1):54-57.

[10]郭軍海. 基于最優(yōu)節(jié)點(diǎn)樣條逼近的觀測(cè)數(shù)據(jù)平滑方法[J]. 中國(guó)空間科學(xué)技術(shù)，2000(3)：43-48.

GUO Jun-hai. A Measurements Smoothing Method Based on Optimal Knots Spline Approximation[J]. Chinese Space Science and Technology, 2000(3):43-48.

[11]潘建平，趙文策，甘友誼.基于測(cè)量數(shù)據(jù)誤差特性求解軌道參數(shù)的方法[J].飛行器測(cè)控學(xué)報(bào)，2009，28(1)：52-55.

PAN Jian-ping, ZHAO Wen-ce, GAN You-yi. Computation of Trajectory Parameters Based on Error Characteristics of Measurement Data[J]. Journal of Spacecraft TT&C Technology, 2009, 28(1):52-55.

[12]趙文策，潘建平，陳偉利.基于彈道動(dòng)力特性考慮的不完全測(cè)量數(shù)據(jù)處理方法[J].飛行器測(cè)控學(xué)報(bào)，2006，25(6)：64-68.

ZHAO Wen-ce, PAN Jian-ping, CHEN Wei-li. The Processing Method of Incomplete Instrumentation Data Based on Consideration of Trajectory Dynamic Characteristics[J]. Journal of Spacecraft TT&C Technology, 2006, 25(6):64-68.

[13]王召剛，袁林，玄志武.PSO節(jié)點(diǎn)尋優(yōu)的樣條逼近微分[J].飛行器測(cè)控學(xué)報(bào)，2012，31(2)：45-48.

WANG Zhao-gang, YUAN Lin, XUAN Zhi-wu. Spline Approximation Diferential with PSO Knot Placement Majorization[J]. Journal of Spacecraft TT&C Technology, 2012, 31(2):45-48.

[14]李鄧化，李金鰲，龐美楓，等.基于慣性因子的混沌粒子群優(yōu)化算法[J].北京信息科技大學(xué)學(xué)報(bào)，2012，27(5)：7-10.

LI Deng-hua, LI Jin-ao, PANG Mei-feng, et al. Research on Chaos Particle Swarm Optimization Algorithm Based on Inertia Weight[J]. Journal of Beijing Information Science and Technology University, 2012, 27(5):7-10.

[15]胥小波，鄭康峰，李丹，等.新的混沌粒子群算法[J].通信學(xué)報(bào)，2012，33(1)：24-30.

XU Xiao-bo, ZHENG Kang-feng, LI Dan, et al. New Chaos Particle Swarm Optimization Algorithm[J]. Journal on Communications, 2012, 33(1):24-30.

Method of Spline Differential Flatness for Position Parameters of Spacecraft with Characteristic Points

CHEN Jing，LI Guo-you，HUANG Wei,ZHANG Xin-lan

(PLA，No. 95948 Troop, Gansu Lanzhou 732750, China)

To obtain the accurate velocity parameters of spacecraft, according to the characteristics of spacecraft orbit with characteristic points, a differential flatness algorithm based on optimal knots spline is proposed. Firstly, the math model of spline differential flatness for position parameters of spacecraft orbit with characteristic points is put forward. Secondly, the recognition method for characteristic points is introduced. Then, the spline differential flatness algorithm based on chaos particle swarm optimization is set up. Finally, the simulation results in typical cases are supplied. The simulation results show that this method can reduce the influence of truncation error and random error, and the accuracy of velocity parameters of spacecraft orbit with characteristic points can be improved.

spacecraft; position parameters; differential flatness; characteristic point; spline knots; particle swarm optimization

2015-07-30；

2016-02-30

陳晶(1978-)，男，遼寧錦州人。工程師，學(xué)士，從事靶場(chǎng)光電測(cè)量總體論證研究。

通信地址：732750蘭州市27支局16信箱1號(hào)E-mail：xdliguoyou@163.com

10.3969/j.issn.1009-086x.2016.04.011

V556.3;TP391.9

A

1009-086X(2016)-04-0061-05