靶場光電外測數據修正及彈道精確復現方法*

魏宗康，郭鎮凈，高榮榮

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

0 引 言

外彈道測量及數據處理是導彈和航天器飛行試驗工程的重要組成部分，它對于保障導彈、航天器試驗的完成和促進其技術發展具有重要的作用[1]。隨著導彈飛行試驗的射程越來越遠，測量精度的要求越來越高，由于光學經緯儀具有不受“黑障區”和地面雜波干擾影響等優點[2]，在導彈、運載火箭飛行試驗的初始段和再入段測量中，光測設備現在仍是主要測量手段之一。但光電經緯儀存在測量距離較近和易受天氣影響等缺點，需要與無線電測量設備相輔相成，共同完成飛行試驗的測量任務。

為克服光電經緯儀測量距離較近的缺點，以獲取更高精度的飛行軌跡，導彈試驗靶場彈道測量設備的種類和數量不斷增加。在雷達外測交接工作前，采用2臺以上光電經緯儀分段測量導彈或者運載火箭的初始段飛行軌跡。但是，不同站點的測量設備在交接過程中引起彈道位置測量參數存在跳臺階的現象。為避免該誤差，交接飛行段采取使用2臺或2臺以上光電經緯儀進行交匯定位的手段，通過不同設備之間的交叉檢驗識別某設備的野值數據，以提高彈道參數的估計精度[3]。

外測數據事后處理的結果主要用于評定導彈和運載火箭的性能和精度、分離制導系統工具誤差等[4]。采用光學或無線電測量設備獲取的飛行器彈道參數一般為位置數據，要想獲得飛行器的速度和加速度，通常采用中心微分平滑方法及相應公式。但是，經過微分后得到的速度噪聲過大，掩蓋了初始段慣性導航的位置誤差和速度誤差。由于外測數據中含有隨機誤差，導致解算結果存在一定程度的波動，因此必須使用平滑濾波算法對解算結果進行處理，才能得到接近真實的飛行彈道參數[5]。常用的平滑濾波主要包含卡爾曼濾波[6]及多項式濾波[7-8]，但由于前者用于實時數據處理，要求精確的運動模型以及部分先驗知識，有時在應用過程中無法滿足該條件；而后者則用于事后處理，基于數據驅動建立參數回歸模型，將數據處理問題轉化為參數估計問題，以提高目標彈道參數精度[9]。常用的濾波算法有中值濾波和低通濾波等，其中應用最廣泛的是卷積平滑算法[10]、樣條插值算法[11]、小波分析[12]和經驗模態分解（Empirical Mode Decomposition，EMD）[13]。但這些估計方法均未討論其在平滑和濾波過程中對重要參數的設定準則，如何科學地設定算法的參數，使其獲得最優的平滑與濾波效果，是一個有待研究的問題[14]。

因此，本文分析了不同站點的光電測量設備在交接過程中引起彈道位置偏差的原因，主要包括站址位置誤差、時間不對齊誤差和隨機誤差三大類，以及噪聲過大帶來的影響。并提出了一種多站交接點位置偏差分段修正的方法，以減小設備交接造成的位置誤差和時間不對齊誤差，再結合多項式擬合的方法降低噪聲對速度的影響，從而得到連續平滑的彈道軌跡，實現彈道軌跡的精確復現。

1 光電經緯儀數據測量模型與誤差產生分析

1.1 光電經緯儀測量模型

在靶場測試飛行器初始段彈道時，采用多站光電經緯儀分段接續測試的方式，各單站分別負責測量不同段的飛行軌跡。光電經緯儀是在光學經緯儀上加裝激光測距系統構成。光學經緯儀為三軸（垂直軸、水平軸、視準軸）地平裝置，水平軸和視準軸可以繞垂直軸在水平面內旋轉，視準軸繞垂直軸旋轉的角度由裝在垂直軸上的碼盤給出（相對某一基準方位），稱為方位角；視準軸繞水平軸旋轉的角度由裝在水平軸上的碼盤給出（水平面為零基準），稱為俯仰角或高低角。

脈沖激光測距通過精確測量激光脈沖信號在測距系統與被測目標之間往返所需的時間來測定距離。由激光發射系統產生并發射峰值功率高、束散角小的激光脈沖，對主波脈沖和回波脈沖之間的時間進行精密測量，并換算成距離值輸出。在主動段測量時，為了提高光電經緯儀的作用距離和測量精度，在導彈和運載火箭上經常安裝激光合作目標[3,15]。

1.1.1 單站測量

圖1 光電經緯儀單點測量示意圖Fig.1 Schematic diagram of photoelectric theodolite single point measurement

采用單個光電經緯儀測量飛行器位置M的示意圖如圖1所示。圖1中，OXYZ為發射坐標系，X軸朝向北，Z軸朝向東，O0點為光電經緯儀站址的位置，M點為導彈或者航天器的位置，則M點在發射坐標系中的投影分量為[15]

（1）

式中，R為光電經緯儀測量的導彈或者航天器相對站址的距離；E和A分別為光電經緯儀測量的高低角和方位角；x0、y0、z0為站址在發射坐標系中的位置分量；x、y、z為由光電經緯儀計算的導彈或航天器在發射坐標系中的位置分量。

1.1.2 多站測量

導彈或運載火箭發射起飛后，由2臺以上光電經緯儀分段測量初始段的飛行軌跡（如圖2所示）。圖2中，采用了3臺光電經緯儀，站址分別為O1（x01，y01，z01），O2（x02，y02，z02），O3（x03，y03，z03），分別負責OA段彈道、AB段彈道和BC段彈道的測量。

圖2 光電經緯儀分段測量彈道示意圖Fig.2 Schematic diagram of segmented ballistic measurement with photoelectric theodolite

當導彈或運載火箭在OABC段彈道飛行時，由這3臺光電經緯儀分段接續進行測量。其中，OA段彈道中M1點、AB段彈道中M2點和BC段彈道中M3點的位置在發射坐標系OXYZ中的投影分量為

（2）

式中，i=1,2,3。

當得到OA段彈道、AB段彈道和BC段彈道的測量數據后，把三段數據按時間序列進行排列，即可得到彈道OABC段的彈道測量值。

1.2 光電經緯儀外測數據誤差產生分析

采用多站測量方式時，會產生以下兩個問題：

1）不同站點的測量設備在交接過程中，彈道位置測量參數存在跳臺階的現象。

采用3臺光電經緯儀測量的導彈初始段彈道如圖3所示，可以看出，測量數據分別在A點、B點交接過程中，彈道發生跳變。

圖3 光電經緯儀分段測量的初始段彈道Fig.3 Initial trajectory measured by photoelectric theodolite

由于3個速度分量vx、vy、vz未發生突變，對三者分別積分后得到一組新的彈道位置數據x′、y′、z′，可得δx=x-x′、δy=y-y′、δz=z-z′，如圖4所示。從圖4中可以看出，δx和δz在發生跳變后其量值近似保持為常值，而δy的值會隨時間單方向漂移。

圖4 外測位置誤差Fig.4 External measurement position error

根據式（2）得到位置誤差方程為

（3）

其中，δ表示差值。

從式（3）可以看出，光電經緯儀外測誤差可以歸為三類：

①常值誤差。主要為光電經緯儀的站址誤差，在A點、B點交接過程中產生的位置誤差主要是由站址O1、O2和O3的定位誤差引起的。該誤差為常值，可以補償。

②時間不對齊誤差。主要是由不同光電經緯儀之間的測量時間不同，以及對位置測量值進行微分時的時間延遲等。該項誤差引起的位置誤差與速度之間成正比，也可進行補償。

③隨機誤差。該項誤差主要是光電經緯儀的測角誤差、激光測距儀的測距誤差和測量數據的量化誤差等，不可補償。

2）速度噪聲過大，掩蓋了慣性導航的速度偏值。

對比圖5中Z軸的z（藍色曲線）和z′（紅色曲線）可以看出，外測位置由于量化誤差引起噪聲偏大，其精度不能如實反映導彈的運動狀態。

圖5 外測位置與速度積分得到的位置對比Fig.5 Comparison between measured position and position obtained by velocity integration

而用外測速度數據時，又面臨著一個新問題，即噪聲水平掩蓋了慣性導航的速度偏值。圖6所示為慣性導航的速度減去外測速度之后的速度差，該速度差的最大幅值接近1m/s。雖然慣性導航的速度誤差會隨時間而積累，但在該時間段的遙外測速度差遠大于慣性導航的速度誤差。因此，需要對外測速度進行平滑濾波，以能夠如實反映慣性導航速度誤差的變化趨勢。

圖6 慣性制導遙外測速度差Fig.6 Velocity difference of inertial guidance remote measurement

2 多站交接點位置偏差修正及外測彈道精確復現方法

為消除光電經緯儀的外測誤差，本文提出了一種多站交接點位置偏差修正方法，旨在通過迭代消除均值差的方式，消除測量設備交接過程中產生的常值誤差和時間不對齊誤差，從而達到修正位置偏差的目的。并且給出了兩種常用的去噪平滑的數據修正方法，以精確復現飛行器的彈道。

2.1 多站交接點位置偏差修正方法

對于由站點位置誤差引起的光電外測跳變誤差，本文針對可以補償的常值位置誤差和時間不對齊誤差，提出了一種多站交接點位置偏差修正方法，通過分段迭代的手段消除均值差，從而修正位置偏差。具體算法流程如下：

（1）以第一個交接站點O2處前后的彈道OAB為例。首先，把OAB段彈道分為三段OA′、A′A、AB。假設該段彈道在A點跳變是由O2的定位誤差引起，位置誤差分量分別為δx02、δy02、δz02，三者為常值。以δx02為例，OAB段位置數據為[XOA′;XA′A;XAB]，位置誤差為[δXOA′;δXA′A;δXAB]，速度數據為[vxOA′;vxA′A;vxAB]。XOA′和δXOA′分別為第一臺光電經緯儀的測量值及其測量誤差；XA′A和δXA′A分別為交接過程的測量值及其測量誤差；XAB和δXAB為第二臺光電經緯儀的測量值及其測量誤差。

（2）其次，根據彈道中速度分量vx的數量級大小，考慮誤差補償類型。當速度數量級較小可以忽略不計時，進入步驟（2.1）；否則進入步驟（2.2）。

（2.1）當該段彈道中速度分量數量級可以忽略不計時，只考慮位置常值誤差。

①首先，給定δx02的初值mx0。對δXAB統一減去mx0，可得數據集[δXOA′;δXA′A;δXAB-mx0]。分別計算δXOA′的均值為mx1，δXAB-mx0的均值為mx2，得到一組新的誤差[δXOA′;δXA′A;δXAB-mx0-mx2+mx1]。其中，δx02=mx0+mx2-mx1。

②其次，給定δXA′A的初值nx0。對δXA′A統一減去nx0，可得數據集[δXOA′;δXA′A-nx0;δXAB-mx0-mx2+mx1]。計算δXA′A-nx0的均值為nx2，得到一組新的誤差[δXOA′;δXA′A-nx0-nx2;δXAB-mx0-mx2+mx1]。

③最后，檢驗誤差[δXOA′;δXA′A-nx0-nx2;δXAB-mx0-mx2+mx1]各段的均值。如果基本相同，得到修正后的外測值為[XOA′;XA′A-nx0-nx2;XAB-mx0-mx2+mx1]。一般情況下，經過上述計算后即可得到滿足要求的結果。如果精度要求較高時，可采用迭代計算。

（2.2）當該段彈道中速度分量數量級不可忽略時，需要同時考慮站點位置誤差和時間不對齊誤差兩部分。

①首先，根據下面的誤差公式，以δXAB和vxAB測量值為例，采用最小二乘算法估計出位置偏差量δrx和時間不對齊誤差量Δt。

δX=δrx+vxΔt

（4）

②其次，把交接設備前的速度數據vxOA′代入式（4），可得

（5）

由此，可求得

（6）

進而得到一組新的外測值[XOA′-vxOA′Δt-δrx;XA′A-vxA′AΔt-δrx;XAB-δx02-vxABΔt-δrx]。

③然后，仍需對δXA′A進行修正。給定δXA′A的初值nx0，對δXA′A統一減去nx0，可得數據集[XOA′-vxOA′Δt-δrx;XA′A-vxA′AΔt-δrx-nx0;XAB-δx02-vxABΔt-δrx]。計算δXOA′的均值為nx1，可得到一組新的誤差為[δXOA′-vxOA′Δt-δrx;δXA′A-vxA′AΔt-δrx-nx0-nx1;δXAB-δx02-vxABΔt-δrx]。

④最后，檢驗誤差[δXOA′-vxOA′Δt-δrx;δXA′A-vxA′AΔt-δrx-nx0-nx1;δXAB-δx02-vxABΔt-δrx]各段的均值。如果基本相同，得到修正后的外測值為[XOA′-vxOA′Δt-δrx;XA′A-vxA′AΔt-δrx-nx0-nx1;XAB-δx02-vxABΔt-δrx]。一般情況下，經上述計算后即可得到滿足要求的結果。如果精度要求較高時，可采用迭代計算。

（3）然后，進入后續彈道ABC段修正過程。該段彈道在B點跳變是由O3的定位誤差引起，位置誤差分量分別為δx03、δy03、δz03，三者為常值。如果速度數量級較小可忽略不計時，借鑒OAB段彈道的步驟（2.1）的位置常值誤差修正方法；否則,需要根據速度修正該段彈道的時間不對齊誤差和常值誤差。流程如下：

①首先，根據步驟（2）得到該段數據誤差為[δXOA′-vxOA′Δt-δrx;δXA′A-vxA′AΔt-δrx-nx0-nx1;δXAB-δx02-vxABΔt-δrx;δXBB′;δXB′C]，把該段數據的速度數據vxBC代入式（4），可得

（7）

進而得到誤差[δXOA′-vxOA′Δt-δrx; δXA′A-vxA′AΔt-δrx-nx0-nx1; δXAB-δx02-vxABΔt-δrx;δXBB′-vxBB′Δt-δrx; δXB′C-vxB′CΔt-δrx]。

②其次，分別計算δXAB-δx02-vxABΔt-δrx的均值為mx1，δXB′C-vxB′CΔt-δrx的均值為mx2，得到一組新的誤差[δXOA′-vxOA′Δt-δrx;δXA′A-vxA′AΔt-δrx-nx0-nx1;δXAB-δx02-vxABΔt-δrx;δXBB′-vxBB′Δt-δrx;δXB′C-vxB′CΔt-δrx-mx2+mx1]。其中，δx03=δX'BC+mx2-mx1。

③最后，檢驗誤差[δXOA′-vxOA′Δt-δrx;δXA′A-vxA′AΔt-δrx-nx0-nx1;δXAB-δx02-vxABΔt-δrx;δXBB′-vxBB′Δt-δrx;δXB′C-vxB′CΔt-δrx-mx2+mx1]各段的均值。如果基本相同，得到修正后的外測值為[XOA′-vxOA′Δt-δrx;XA′A-vxA′AΔt-δrx-nx0-nx1;XAB-δx02-vxABΔt-δrx;XBB′-vxBB′Δt-δrx;XB′C-vxB′CΔt-δrx-mx2+mx1]。一般情況下，經上述計算后即可得到滿足要求的結果。如果精度要求較高時，可采用迭代計算。

2.2 測量噪聲平滑方法

光電經緯儀數據處理時，常利用由位置參數微分中心平滑的方法及相應公式，獲取導彈的飛行速度和加速度參數。對于非關機點附近的主動段彈道、自由段彈道和再入段彈道，通常使用速度二階中心平滑和加速度三階中心平滑公式處理光電經緯儀的數據。處理主動段彈道數據時,平滑區間常取為1s或2s;處理自由段彈道數據時,平滑區間一般不超過20s;再入段彈道的平滑區間為1～2s。對于關機點附近的主動段彈道，由于彈道變化劇烈，由位置參數平滑求速時，采用四階中心平滑公式，平滑區間同二階中心平滑[4,9,16]。

但是微分得到的數據存在跳變和噪聲過大的現象，掩蓋了慣性導航的速度偏值，因此需要對外測速度進行平滑濾波，以能夠如實反映慣性導航速度誤差的變化趨勢。下面給出兩種常用的濾波方法。

2.2.1 低通濾波方法

低通濾波器作為一種廣泛應用的濾波器，具備保留低頻信號、抑制高頻信號的特點，并且性能穩定，易于調節，工程上常用于信號處理、數據傳輸和抑制干擾。利用這一特點，可使用在外測彈道速度和加速度參數處理中，由于彈道速度是由位置微分后得到的，周期為2s。采用如下的二階低通濾波器

（8）

（a） f=1Hz

（b） f=0.1Hz圖7 低通濾波效果對比Fig.7 Comparison of low-pass filtering effect

對速度進行再濾波時，如果轉折頻率大于1Hz，則濾波無效果；而如果轉折頻率過小，則響應能力變弱，曲線失真。圖7所示為分別采用f=1Hz和f=0.1Hz對vz進行濾波的結果，圖中藍色曲線為vz濾波前的值，紅色曲線為vz濾波后的值。可以看出，f=1Hz時，噪聲幅值無變化；而f=0.1Hz時，噪聲雖然衰減，但在信號快速變化時響應能力不夠。因此，在后續分析中不選用濾波法，而采用曲線擬合方法。

2.2.2 多項式擬合算法

多項式擬合算法也叫Savitzky-Golag（SG）平滑算法，是由Savitzky和Golag提出的基于最小二乘原理的多項式平滑算法[8,10,17]。該算法的關鍵在于矩陣算子的求解。假定平滑濾波滑動窗口長度（數據點的個數）為k，采用n次多項式對窗口內的數據進行描述[14]，如式（9）所示。通常情況下，滑動窗口的長度和多項式擬合階次的選擇來源于工程經驗。

（9）

式中，i=1,2,…,k。

將k個測量點依次代入多項式方程中，得到了k個方程組，將其轉化為矩陣形式為

（10）

對以上矩陣方程采用最小二乘估計，可得到擬合參數ai（i=1,2,…,k）的估值。

3 試驗結果及分析

3.1 多站交接點位置偏差修正結果

針對圖4所示的由站點位置誤差引起的光電外測跳變誤差，采用本文方法進行修正：

（1）OAB段彈道修正

該段彈道在A點跳變是由O2的定位誤差引起，位置誤差分量分別為δx02、δy02、δz02，三者為常值。由于Y方向的速度分量vy相對vx和vz的值大1～2個數量級，需單獨修正彈道誤差。對Y方向彈道采用最小二乘法估計出站址位置誤差δry=-1.266m和時間不對齊誤差Δt=7.933×10-3s。修正后OAB段外測位置誤差如圖8所示，對比圖4可以看出，這種修正方法既可辨識出站址常值誤差，也可辨識出時間不對齊誤差，經過補償后實現了對OAB段彈道的精確復現。

圖8 修正后的外測位置誤差Fig.8 Corrected external measurement position error

（2）OABC段彈道修正

圖9 修正后的外測位置誤差Fig.9 External measurement position error after correction

該段彈道在B點跳變是由O3的定位誤差引起，位置誤差分量分別為δx03、δy03、δz03，三者為常值。針對圖8所示的由站點位置誤差引起的光電外測跳變誤差，修正方法借鑒OAB段彈道的修正方法，但比較特殊的是δy03值的精確估計及彈道修正，仍需要修正時間不對齊引起的誤差。修正后OABC段外測位置誤差如圖9所示，對比圖4和圖8可以看出，修正后的測試數據不再發生跳變，經過站址誤差標定和補償后實現了對OABC段彈道的精確復現。

經過多站交接點位置偏差修正方法，實現了站點交接引起的彈道位置偏差和時間不對齊偏差。修正后的彈道如圖10中紅色曲線所示，相對于修正前的藍色曲線，彈道更加連續。

圖10 修正前后的初始段彈道Fig.10 Initial trajectory before and after correction

3.2 多項式擬合平滑噪聲結果

本文滑動窗口長度取光電測試數據的長度，取k=1120，n=4，對速度測量值vx、vy、vz進行平滑擬合，計算得到3個方向速度的擬合參數分別為

結果如圖11所示，其中，紅色曲線為速度測量值;藍色曲線為擬合結果；黑色曲線為擬合殘差。

利用上述新的外測速度值，可得到一組新的遙外測速度誤差，如圖12中紅色曲線所示，相比修正前的藍色曲線可以看出，數據更加平滑。

最終，對速度進行積分后可得位置信息，如圖13中綠色曲線所示。可以看出，通過修正位置偏差和多項式擬合的方法,消除噪聲后得到的彈道更加平滑準確，既保證了彈道曲線的連續性，又可滿足彈道的平滑性。

圖11 速度擬合值及擬合殘差Fig.11 Velocity fitting value and fitting residual

圖12 慣性制導遙外測速度差Fig.12 Velocity difference of inertial guidance remote measurement

圖13 修正偏差和平滑噪聲前后的初始段彈道Fig.13 Initial trajectory before and after correcting bias and smoothing noise

4 結 論

在導彈試驗靶場采用多套光電外測設備測量外彈道參數時，不同設備在交接工作的過程中會引起位置偏差，后續利用外測位置數據獲得速度和加速度時，由于微分會存在噪聲過大的現象，從而導致測量彈道不準確，不利于后續的試驗分析。因此，本文提出了一種多站交接點位置偏差的修正方法，旨在通過迭代消除均值差的方式，消除測量設備交接過程中產生的常值誤差和時間不對齊誤差，從而達到修正位置偏差的目的，使得彈道曲線更加連續。同時，利用多項式擬合的方法緩解由微分計算速度值而帶來的噪聲過大的問題，使得彈道曲線更加平滑，從而較好地復現彈道飛行軌跡。但是，濾波過程中如何從理論上科學地設定算法的重要參數，使其獲得最優的平滑效果，以及如何應用平滑后的運動參數，實現導彈的精度評定、制導系統工具誤差分離等工作，需要進一步的研究和驗證。