基于樣條函數表征目標運動軌跡事后數據融合方法研究

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

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

0 引言

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

1 數據融合模型

1.1 聯合誤差方程建立與解算

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

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

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

定位測元誤差方程

距離測元誤差方程

方位角測元誤差方程

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

徑向速度測元誤差方程

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

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

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

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

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

1.2 軌跡誤差理論分析

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

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

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

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

1.3 樣條節點優化擬合

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

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

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

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

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

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

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

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

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

2 應用與分析

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

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

結合以上聯合誤差方程建立和最優節點搜索方法進行靶機飛行軌跡數據融合計算。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結論

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

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

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

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

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

References)

[1]朱武宣，高耀文.基于樣條約束“EMBET”的再入軌道測量數據融合方法[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]王正明，易東云，周海銀.彈道跟蹤數據的校準與評估[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 和陸基測量的外彈道數據事后融合方法研究[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]宮志華，董立濤，徐旭.外彈道測量數據融合方法研究[J].兵器試驗，2012，261(3):1 －6.GONG Zhi－hua，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].北京:電子工業出版社，2007.WANG Zheng－lin，GONG Chun，HE Qian.Master matlab scientitic calculation[M].Beijing:Publishing House of Electronics Industry，2007.(in Chinese)

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

[8]周海銀.空際目標跟蹤數據的融合理論和模型研究及應用[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)