宮志華，周海銀，郭文勝，徐旭

(1.中國白城兵器試驗中心，吉林 白城137001;2.國防科學技術大學 理學院，湖南 長沙410073)

0 引言

針對常規(guī)武器系統(tǒng)發(fā)射目標機動性強、尺寸小、飛行高度低的特點，迫使測控系統(tǒng)向著一個多設備、多體制(光學、雷達、遙感)、多冗余的綜合測控體系發(fā)展，也驅動著數(shù)據(jù)處理性能精準化程度要求日益提高。傳統(tǒng)、獨立、互不交融的數(shù)據(jù)處理方法不再滿足數(shù)據(jù)處理要求，在此情形下，如何進行多冗余多結構異類測元數(shù)據(jù)綜合處理，以獲得更可靠的試驗結果一直是學術界一個研究的熱點。傳統(tǒng)運動目標軌跡參數(shù)融合解算方法認為，運動目標在時序上不具有任何關聯(lián)性，即在每一個采樣時刻進行目標軌跡參數(shù)融合解算，這種方法沒有充分利用測量數(shù)據(jù)之間的時序關聯(lián)性[1]，使設計矩陣龐大，耗費運算資源，效率低。針對這一問題，本文結合機動、小型運動目標特點，提出采用基于自由節(jié)點樣條函數(shù)來描述目標運動軌跡并轉換，實現(xiàn)以其為基函數(shù)對測控設備測元進行表征，將冗余多源異類測元數(shù)據(jù)、多采樣時刻的測量方程聯(lián)合求解目標軌跡參數(shù)的融合方法。該方法充分考慮運動目標全程軌跡特征點和時序關聯(lián)性，極大減小由數(shù)學模型表征目標軌跡造成的截斷誤差，大量壓縮待估參數(shù)數(shù)量，提高目標軌跡參數(shù)估計精度和計算穩(wěn)定度，同時對各測元潛在的系統(tǒng)誤差進行多循環(huán)探查和估計，實現(xiàn)系統(tǒng)誤差自校準的目的。

1 數(shù)據(jù)融合模型

1.1 聯(lián)合誤差方程建立與解算

考慮不同測控設備測量原理，可以提供的獨立測元種類有多種，這里只列舉常用的3 種:GPS、GLONASS、北斗空間定位測元(局部坐標(x，z，y)，由原始大地坐標(緯度b，經度l，大地高h)轉換得到)、光學設備測元(方位角α，俯仰角θ)、脈沖雷達測元(距離r，方位角α，俯仰角θ)和連續(xù)波雷達測元(距離r，方位角α，俯仰角θ，徑向速度v).下面，就以靶場常見的7 種獨立測元(x，z，y，r，α，θ，v)為例，具體分析基于樣條函數(shù)表征測元的聯(lián)合誤差方程的建立和解算過程。

對(1)式進行一階微分，可得到(2)式，

關于各測元誤差方程的建立，除定位測元(x，z，y)外，其他測元(r，α，θ，v)均為非線性方程，需按泰勒級數(shù)展開轉化為線性方程。則以(1)式、(2)式為基函數(shù)表征各測元誤差方程公式推導總結如下:

定位測元誤差方程

距離測元誤差方程

方位角測元誤差方程

式中:αT為判象限角。俯仰角測元誤差方程

徑向速度測元誤差方程

(3)式～(7)式中:(xi，zi，yi，ri，αi，θi，vi)為各測元實測數(shù)據(jù);(sx，sz，sy，sr，sα，sθ，sv)為各測元系統(tǒng)誤差模型，如常值、線性或非線性函數(shù)模型等;(初 值;(x0，z0，y0)為 站 址 已 知 坐 標;ri0=

將(3)式～(7)式各測元誤差方程聯(lián)立，組成目標軌跡聯(lián)合誤差方程，寫成矩陣形式，如(8)式所示。

式中:ζ=(exi，ezi，eyi，eri，eαi，eθi，evi)T為測元殘差向量;H 為與樣條函數(shù)有關的系數(shù)矩陣(亦稱為設計矩陣);X = (βx1，…，βx(mx+n+1)，βz1，…，βz(mz+n+1)，βy1，…，βy(my+n+1))T為與樣條函數(shù)有關的待估參數(shù)向量;G 為系統(tǒng)誤差模型系數(shù)矩陣;C 為系統(tǒng)誤差模型待估參數(shù)向量;η=(xi，zi，yi，ri－ri0+riΔ，αi－αi0+αiΔ，θi－θi0+θiΔ，vi)T為常數(shù)向量。

在參數(shù)估計時，在大數(shù)據(jù)量處理的情況下，相關性和非高斯性對估計結果影響不大[3]，因此采用白噪聲條件下的參數(shù)估計方法，在已知最優(yōu)節(jié)點分布的前提下，由最小二乘法原理解算可以得到待估參數(shù)(樣條擬合函數(shù)系數(shù)和系統(tǒng)誤差模型系數(shù))和其協(xié)方差矩陣解，如(9)式所示，

式中:P 為權值矩陣。在實際求解中，由于初始目標軌跡參數(shù)的近似性和非線性函數(shù)級數(shù)展開帶來的截斷誤差，需要進行疊代計算。將最終解算得到的樣條函數(shù)擬合系數(shù)和系統(tǒng)誤差函數(shù)系數(shù)代入目標軌跡參數(shù)表達式和誤差函數(shù)表達式中，即可得到融合計算后的目標軌跡坐標值和系統(tǒng)誤差值。

1.2 軌跡誤差理論分析

現(xiàn)在考慮融合后目標軌跡坐標的精度估算問題。首先，測元值可寫成(10)式的矩陣形式，

即測元值Y 是由以樣條擬合系數(shù)β 和系統(tǒng)誤差模型系數(shù)φ 為參數(shù)的F 函數(shù)表征，ξ 為測元殘差向量。這樣，將函數(shù)F 在樣條系數(shù)估值處展開[4]，則近似地有

式中:H(t)=(Bx(t)T，Bz(t)T，By(t)T)T為樣條基系數(shù)矩陣;X(t)為目標軌跡坐標的真值。目標軌跡坐標協(xié)方差為

式中:Λ 為由測元殘差ξ 的方差組成的對角陣。

1.3 樣條節(jié)點優(yōu)化擬合

以上采用樣條函數(shù)表征目標軌跡，為減小截斷誤差，必須優(yōu)選樣條節(jié)點分布。關于樣條節(jié)點的選取，首先的一個想法是可以采用等距節(jié)點。等距節(jié)點的實質是相當于采用一個恒定頻率帶寬的濾波器對數(shù)據(jù)進行濾波，但是由于非平穩(wěn)數(shù)據(jù)所含噪聲頻率成分較多，且存在特征點現(xiàn)象，顯然，采用恒定濾波器對數(shù)據(jù)濾波，要達到更好的精度，必然要使用更密的節(jié)點數(shù)以減小濾波器帶寬，這對于節(jié)省待估參數(shù)進行融合解算是不利的[5]，因此，針對機動、非平穩(wěn)目標軌跡特性，必須采用自由節(jié)點樣條函數(shù)來對其表征，并在一定約束條件下對節(jié)點分布進行優(yōu)化搜索。

如(1)式所示，用樣條函數(shù)表征目標軌跡，以x方向為例，可以表示成矩陣形式，如(15)式，

式中:a= (a0，a1，a2，…，amx)T為樣條系數(shù)參數(shù)向量;Y=(x1，x2，…，xn)T為目標軌跡向量;ε =(ε1，ε2，…，εn)T為隨機誤差向量。

上式參數(shù)估計可歸為如下非線性優(yōu)化問題:即求樣條擬合系數(shù)a 和節(jié)點分布TM(M 為節(jié)點個數(shù))，使得

式中:Q(a，TM)=‖Y －X(TM)a‖2=(Y －X(TM)·a)T(Y－X(TM)a).

為定量給出低頻信號與隨機誤差的分頻界限，定義一個BIC 值進行判斷[2]，如(17)式所示，

對于(16)式的多參數(shù)非線性目標函數(shù)，可以采用多種優(yōu)化算法，文獻[5]對這些算法的原理和實現(xiàn)都有詳細的說明。本文采用的是模擬退火法。

需要特別注意的是，特征點必須作為已知樣條節(jié)點使用，參與搜索其他最優(yōu)節(jié)點，以避免產生較大的截斷誤差。特征點的獲得可以采用遙測信息、光學圖像判讀、雷達回波譜分析和目標軌跡加速度曲線[7]分析等多種手段。

圖1 目標軌跡三維坐標樣條函數(shù)擬合BIC 值曲線Fig.1 The BIC values of three－dimensional trajectory based on spline function representation

2 應用與分析

對以上數(shù)據(jù)融合方法的合理性、正確性和實用性進行驗證分析，結合某靶機飛行試驗展開，分別有光學、雷達、遙感等測控設備參試跟蹤靶機，包括紅外光學經緯儀、脈沖雷達和連續(xù)波雷達，其中遙測設備記錄GPS 機載設備星歷原始數(shù)據(jù)，事后提供經過GPS 載波相位差分修正后的高精度靶機軌跡坐標數(shù)據(jù)作為比對真值數(shù)據(jù)，定位精度優(yōu)于0.5 m.

選取8 個獨立原始測元參與靶機軌跡融合計算，包括紅外光學經緯儀兩站測量的4 個獨立測元數(shù)據(jù)(方位角和俯仰角)、脈沖雷達測量的3 個獨立測元數(shù)據(jù)(距離、方位角和俯仰角)和連續(xù)波雷達測量的1 個獨立測元數(shù)據(jù)(徑向速度)。各測元原始數(shù)據(jù)采樣頻率不一致。

結合以上聯(lián)合誤差方程建立和最優(yōu)節(jié)點搜索方法進行靶機飛行軌跡數(shù)據(jù)融合計算。

首先，對8 個獨立測元實測數(shù)據(jù)分別進行已知系統(tǒng)誤差修正、跨零跳點修正、合理性檢驗、數(shù)據(jù)插值和數(shù)據(jù)平滑等預處理項工作，本方法不需要進行地球曲率修正。其中，最后比較重要的一項內容是對各測元實測數(shù)據(jù)進行隨機誤差統(tǒng)計，可參考文獻[6]相關方法。在數(shù)據(jù)融合處理中，各測元權重是不一樣的，隨機誤差小表示精度高，權值相對大;反之，權值相對小。因此，以獨立測元等方差和各測元不相關為依據(jù)，由各測元統(tǒng)計方差組成權值系數(shù)對角矩陣。

其次，對樣條函數(shù)最優(yōu)節(jié)點進行搜索確定，選取目標軌跡初值，可以為光學交會數(shù)據(jù)、雷達定位數(shù)據(jù)或理論軌跡數(shù)據(jù)等。考慮本實例中目標軌跡變化比較平穩(wěn)，因此，先以脈沖雷達定位數(shù)據(jù)為初值，采用等距節(jié)點進行聯(lián)合方程的解算，測元系統(tǒng)誤差不考慮，再用融合后的目標軌跡數(shù)據(jù)進行最優(yōu)節(jié)點的搜索。其中，公共段數(shù)據(jù)采用預處理后的50 ～100 s數(shù)據(jù)，數(shù)據(jù)采樣頻率為10 Hz，則按5 s 等間隔取內節(jié)點數(shù)量為9 個。

然后，分別對初次融合后獲得的目標軌跡進行三次B 樣條函數(shù)擬合，按以上方法搜索樣條函數(shù)最優(yōu)節(jié)點分布，得到目標軌跡三維坐標樣條函數(shù)擬合效果的BIC 判定值，如圖1 所示。

從圖1 可以判斷目標軌跡在三維方向上由樣條函數(shù)擬合的最優(yōu)節(jié)點數(shù)不是9 個，而是10 個，接著將搜索獲得的最優(yōu)節(jié)點分布分別代入(1)式，并依據(jù)以上建立的聯(lián)合誤差方程求解待估參數(shù)。

關于各測元系統(tǒng)誤差模型的確定，在實際解算過程中可以采用以下過程:1)通過理論分析、比對等手段可以先確定一些測元系統(tǒng)誤差的先驗信息;2)無先驗信息情況下，在初次融合解算時，先不考慮系統(tǒng)誤差，查看融合解算后的各測元殘差曲線特征;3)對殘差有明顯趨勢項的測元，依殘差趨勢項特征建模，并重新融合解算，再查看各測元殘差曲線特征，循環(huán)往復，直到所有殘差均值為0;4)各測元系統(tǒng)誤差模型需盡量簡單。依據(jù)的原則是，如果各測元不存在系統(tǒng)誤差或系統(tǒng)誤差建模正確，則融合后各測元對應的殘差應沒有明顯的趨勢項，接近純隨機誤差[5]。

本實例在無誤差融合解算后發(fā)現(xiàn)多數(shù)測元殘差曲線的均值不為0，如圖2 所示，統(tǒng)計結果如表1 所示，表1 中α1、θ1為1 號經緯儀測元，α2、θ2為2 號經緯儀測元，α3、θ3、r3為脈沖雷達測元，v 為連續(xù)波雷達測元。這表明參與數(shù)據(jù)融合的8 個獨立測元中，某些測元必定含有系統(tǒng)誤差，并感染到了其他測元。

圖2 不考慮系統(tǒng)誤差情況下融合解算后各測元殘差曲線Fig.2 The residual of each measuring element based on data fusion without system error

表1 不考慮系統(tǒng)誤差融合解算后各測元殘差統(tǒng)計值Tab.1 The residual statistic of each measuring element based on data fusion without system error

經過不斷地對各測元系統(tǒng)誤差模型的假設建立，并查看融合計算后殘差曲線的均值情況，最后判定兩站光學經緯儀的俯仰角測元和脈沖雷達的3 個測元均含有常值系統(tǒng)誤差。這樣，在最終聯(lián)合方程中需要求解的待估參數(shù)總數(shù)量是47 個，其中，樣條函數(shù)待估系數(shù)參數(shù)是42 個，系統(tǒng)誤差待估參數(shù)是5 個。經過最后的融合解算，得到各測元殘差曲線，如圖3 和表2 所示。從表2 中可以明顯看到，各測元殘差均值都為0;測元的系統(tǒng)誤差值也被探查出來，可用于修正原始測元數(shù)據(jù)。目標軌跡坐標的理論誤差估計值如圖4 所示，統(tǒng)計結果如表3 所示。

圖3 考慮系統(tǒng)誤差情況下融合解算后的各測元殘差曲線Fig.3 The residual of each measuring element based on data fusion with system error

表2 考慮系統(tǒng)誤差融合解算后各測元殘差統(tǒng)計值Tab.2 The residual statistic of each measuring element based on data fusion with system error

表3 融合目標軌跡理論誤差統(tǒng)計結果Tab.3 The theoretical error statistic based on data fusion

下面，為能直觀地表現(xiàn)融合數(shù)據(jù)的質量和效果，以GPS 定位數(shù)據(jù)為真值，與兩站光學經緯儀交會的目標軌跡坐標數(shù)據(jù)和由數(shù)據(jù)融合得到的目標軌跡坐標數(shù)據(jù)進行比對，查看兩組比對誤差，如圖5 和圖6所示，統(tǒng)計結果如表4 所示。

圖5 光學交會目標軌跡數(shù)據(jù)與GPS 定位數(shù)據(jù)比對誤差曲線Fig.5 The error comparison between optical－theodolite and GPS data

圖6 數(shù)據(jù)融合目標軌跡數(shù)據(jù)與GPS 定位數(shù)據(jù)比對誤差曲線Fig.6 The error comparison between fused trajectory data and GPS data

表4 兩種目標軌跡與GPS 數(shù)據(jù)的比對殘差統(tǒng)計結果Tab.4 The residual statistic comparison of two larget trajectories with GPS data

從表4 比對誤差統(tǒng)計結果來看:1)數(shù)據(jù)融合比光學經緯儀交會的目標軌跡坐標數(shù)據(jù)方差和均值都小，更接近真值;2)兩組比對誤差均值都不為0，因為都含有跟蹤部位不一致和大氣折射誤差修正殘差所帶來的誤差，但光學交會數(shù)據(jù)還存在由測元系統(tǒng)誤差造成的交會誤差。

3 結論

基于以上分析，可以得到如下結論:

1)結合運動目標軌跡時序關聯(lián)特性，可以實現(xiàn)以樣條函數(shù)為基函數(shù)對測控設備各種測元進行表征;結合目標運動特征點情況和建立合理目標函數(shù)對樣條函數(shù)最優(yōu)節(jié)點進行優(yōu)化搜索解算，使節(jié)點分布自適應目標運動特性，極大減小截斷誤差。這兩種手段使得聯(lián)合誤差方程的待估參數(shù)被極大壓縮，冗余程度極大增強，計算穩(wěn)定性和效率也提高。

2)建立基于樣條函數(shù)表征目標軌跡的稀疏參數(shù)聯(lián)合誤差方程，在實際計算過程中，對初值準確度要求不高，迭代計算收斂快，解算精度高，這對于需要多次循環(huán)反復進行的融合處理過程將極大地提高處理效率，并對各測元潛在系統(tǒng)誤差具有較敏感的探查能力。因此，該方法實際應用效果明顯。

3)上述針對運動目標軌跡數(shù)據(jù)融合計算方法得到的結果不是最優(yōu)。因為，由于樣條函數(shù)最優(yōu)節(jié)點和權值都需要在融合前確定，由目標軌跡初值或半程融合得到的目標軌跡來約束搜索獲得的節(jié)點分布不是最優(yōu);由高斯—馬爾可夫定理又知，線性方程融合最優(yōu)權值可由統(tǒng)計方差獲得，但非線性方程并不是最優(yōu)[7]。

關于目標運動軌跡數(shù)據(jù)融合研究領域，當前，有專家學者提出在融合模型中將節(jié)點和權值都當作待估參數(shù)，通過融合模型迭代計算改進節(jié)點和權值估計[8]，這種方法由于解空間維數(shù)高，計算量大、非線性程度高和效率低等問題，其工程實用性還有待進一步完善。因此，這一領域的工程研究今后還會不斷深入地進行下去。

References)

[1]朱武宣，高耀文.基于樣條約束“EMBET”的再入軌道測量數(shù)據(jù)融合方法[J].飛行器測控學報，2005，24(6):49 －53.ZHU Wu－xuan，GAO Yao－wen.A data fusion method for re－entry ballistic measurementwith spline restraint EMBET[J].Journal of Spacecraft TT＆C Technology，2005，24(6):49 － 53.(in Chinese)

[2]王正明，易東云，周海銀.彈道跟蹤數(shù)據(jù)的校準與評估[M].長沙:國防科技大學出版社，1999.WANG Zheng－ming，YI Dong－yun，ZHOU Hai－yin.Tracking trajectory data calibrating and evaluating[M].Changsha:National University of Defense Technology Publishing House，1999.(in Chinese)

[3]原野.基于GNSS 和陸基測量的外彈道數(shù)據(jù)事后融合方法研究[D].長沙:國防科學技術大學，2008.YUAN Ye.Research on the post－flight method of exterior ballistic measurement with GNSS data and ground－based TT ＆C network data[D].Changsha:National University of Defense Technology，2008.(in Chinese)

[4]童麗，易東云.非線性融合模型的彈道估計精度評定[J].彈道學報，2002，14(4):1 －5.TONG Li，YI Dong－yun.Assess accuracy of orbit parameters estimated from nonlinear fusion model[J].Journal of Ballistics，2002，14(4):1 －5.(in Chinese)

[5]宮志華，董立濤，徐旭.外彈道測量數(shù)據(jù)融合方法研究[J].兵器試驗，2012，261(3):1 －6.GONG Zhi－h(huán)ua，DONG Li－tao，XU Xu.The method of exterior measured trajectory data fusion[J].Ordance Test，2012，261(3):1 －6.(in Chinese)

[6]王正林，龔純，何倩.精通MATLAB 科學計算[M].北京:電子工業(yè)出版社，2007.WANG Zheng－lin，GONG Chun，HE Qian.Master matlab scientitic calculation[M].Beijing:Publishing House of Electronics Industry，2007.(in Chinese)

[7]朱炬波.不完全測量數(shù)據(jù)建模與應用[D].長沙:國防科學技術大學，2004.ZHU Ju－bo.The incompletely measured data modeling and appling[D].Changsha:National University of Defense Technology，2004.(in Chinese)

[8]周海銀.空際目標跟蹤數(shù)據(jù)的融合理論和模型研究及應用[D].長沙:國防科學技術大學，2004.ZHOU Hai－yin.Researches on theories and models of spatial targets tracking data fusion with applications[D].Changsha:National University of Defense Technology，2004.(in Chinese)