空間碎片激光測距的距離預(yù)報(bào)偏差實(shí)時(shí)修正算法

高健，梁智鵬，韓興偉，2，董雪，安寧，溫冠宇，劉承志，2

（1 中國科學(xué)院國家天文臺(tái)長春人造衛(wèi)星觀測站，長春 130117）

（2 中國科學(xué)院空間目標(biāo)與碎片觀測重點(diǎn)實(shí)驗(yàn)室，南京 210008）

0 引言

衛(wèi)星激光測距技術(shù)（Satellite Laser Ranging，SLR）作為單點(diǎn)精度最高的空間大地測量技術(shù)之一，其測距精度已接近毫米級(jí)，將SLR 技術(shù)應(yīng)用在空間碎片的精密測定軌中能夠使其軌道精度得到大幅提升［3-6］。自2002 年澳大利亞EOS 公司GREENE B 等［7］首次報(bào)告實(shí)現(xiàn)了空間碎片的漫反射激光測距后，國內(nèi)外多個(gè)激光測距臺(tái)站也相繼開展了空間碎片測距實(shí)驗(yàn)，并成功獲取到空間碎片的激光觀測數(shù)據(jù)［8-11］。

空間碎片測距相比合作目標(biāo)測距的主要難點(diǎn)在于極低的漫反射回波接收信噪比和較大的初軌預(yù)報(bào)偏差。近年來針對空間碎片測距的方法創(chuàng)新和技術(shù)改進(jìn)主要集中于提高測距系統(tǒng)的微弱回波探測能力，如近紅外波長測距［12］、多脈沖測距［13-14］、單站發(fā)射多站接收［15-17］、超導(dǎo)納米線單光子探測［18-19］等多種技術(shù)手段的組合與疊加。然而測距成功率由探測概率和虛警概率共同決定，為了抑制天空背景噪聲，使虛警概率保持在較低水平，普遍采用時(shí)間濾波的距離門控技術(shù)，其依賴于高精度的距離預(yù)報(bào)值來計(jì)算精確的探測器開門時(shí)刻，而較大的距離預(yù)報(bào)偏差將加大回波捕獲難度，降低回波搜索效率，甚至無法在目標(biāo)過境弧段內(nèi)完成測距過程。

在提升空間碎片的初軌預(yù)報(bào)精度以及預(yù)報(bào)偏差的實(shí)時(shí)改進(jìn)方面，通過利用已產(chǎn)生的短弧長稀疏激光數(shù)據(jù)或融合光學(xué)測角數(shù)據(jù)能夠提升有限數(shù)量目標(biāo)的測站預(yù)報(bào)精度［20-21］，但“后處理”的改進(jìn)方式以預(yù)先獲取激光數(shù)據(jù)為前提，無法滿足激光測距對于高精度預(yù)報(bào)的高實(shí)時(shí)性需求。通過解算并應(yīng)用沿跡方向的時(shí)間偏差能夠在一定程度上間接減小距離偏差［22］，奧地利Graz 觀測站在空間碎片白天測距中即應(yīng)用時(shí)間偏差實(shí)時(shí)解算技術(shù)來減小距離搜索范圍［23］，但仍無法對視線方向的距離預(yù)報(bào)偏差直接進(jìn)行修正。

本文借鑒衛(wèi)星編隊(duì)飛行［24］的概念，將空間碎片的預(yù)報(bào)軌道與真實(shí)軌道看作鄰近軌道，從空間目標(biāo)軌道相對運(yùn)動(dòng)的角度來解決TLE 預(yù)報(bào)偏差的改進(jìn)問題。以光學(xué)測角信息作為觀測量，通過預(yù)報(bào)軌道上虛擬預(yù)報(bào)目標(biāo)與真實(shí)目標(biāo)的視位置角度偏差來反演距離預(yù)報(bào)偏差的演化過程。以卡爾曼濾波為狀態(tài)估計(jì)手段，構(gòu)建了近圓軌道目標(biāo)的距離預(yù)報(bào)偏差實(shí)時(shí)修正算法，并通過數(shù)值仿真對算法的有效性進(jìn)行了驗(yàn)證。

1 理論模型建立

空間碎片預(yù)報(bào)目標(biāo)與真實(shí)目標(biāo)分別沿各自軌道繞地心運(yùn)行，其相互間的軌道運(yùn)動(dòng)狀態(tài)服從相對運(yùn)動(dòng)動(dòng)力學(xué)規(guī)律?？臻g碎片相對于觀測站的距離預(yù)報(bào)偏差作為預(yù)報(bào)偏差的視線方向分量，無法直接通過光學(xué)測角手段進(jìn)行測量和修正，而處在鄰近軌道的空間碎片預(yù)報(bào)目標(biāo)與真實(shí)目標(biāo)的視位置角度偏差則可通過光學(xué)測角手段實(shí)時(shí)獲取，且相對運(yùn)動(dòng)動(dòng)力學(xué)方程為預(yù)報(bào)偏差的角度分量與距離分量之間的相互轉(zhuǎn)化提供了數(shù)學(xué)基礎(chǔ)。因此，以相對運(yùn)動(dòng)動(dòng)力學(xué)方程為理論建模依據(jù)，通過建立相應(yīng)的空間幾何模型、狀態(tài)演化模型與角度觀測模型，即可實(shí)現(xiàn)距離偏差真值的求解。

1.1 空間幾何模型

將空間碎片預(yù)報(bào)目標(biāo)作為主星，真實(shí)目標(biāo)作為從星，構(gòu)建圖1 所示的空間幾何模型，其包含了主星、從星、地心、觀測站之間的空間位置幾何關(guān)系。圖1 中，O點(diǎn)為預(yù)報(bào)目標(biāo)質(zhì)心，P點(diǎn)為真實(shí)目標(biāo)質(zhì)心，E點(diǎn)為地心，S點(diǎn)為觀測站，Oxyz為固連在主星O上的星基軌道坐標(biāo)系，x軸和y軸在軌道平面內(nèi)，x軸方向由地心E指向主星O，y軸垂直于x軸并指向主星O的運(yùn)動(dòng)方向，z軸指向軌道平面法向，與x軸和y軸構(gòu)成右手直角坐標(biāo)系。γ為由觀測站S指向主星O的單位矢量，同時(shí)定義γ1為單位矢量γ在軌道平面內(nèi)的投影與x軸的夾角，γ2為單位矢量γ與z軸的夾角。

(1)廬山交通發(fā)展快速。2017年7月，廬山索道正式開通，并在開通后的首月共接待游客超過20萬人次，受到了廣大游客的喜愛。廬山索道項(xiàng)目也是九江市旅游的重點(diǎn)工程，僅8分鐘就可以上下廬山，它的建造得到了江西省委、省政府、市委、市政府和廬山管理局的高度重視和支持，有效緩解了廬山旅游旺季的交通壓力，為游客上廬山提供了極大的便利。

1.2 狀態(tài)演化模型

在相對運(yùn)動(dòng)動(dòng)力學(xué)方程中，Clohessy-Wiltshire（C-W）方程［25］是描述近圓軌道目標(biāo)相對運(yùn)動(dòng)的簡潔形式?？紤]距離偏差修正過程僅針對空間碎片的單次短時(shí)過境，其攝動(dòng)項(xiàng)影響遠(yuǎn)小于測量誤差，因此，將C-W方程作為求解預(yù)報(bào)偏差矢量演化過程的數(shù)學(xué)模型，其解析表達(dá)式如式（1），其中(x，y，z)為真實(shí)目標(biāo)在預(yù)報(bào)目標(biāo)質(zhì)心坐標(biāo)系下的相對位置空間坐標(biāo)，對應(yīng)于預(yù)報(bào)偏差矢量的空間狀態(tài)，n為其繞地心運(yùn)行的平角速度。

由式（2）和式（3）的C-W 方程解的形式可知，預(yù)報(bào)目標(biāo)與真實(shí)目標(biāo)在任意時(shí)刻t的相對位置和相對速度狀態(tài)僅由相對運(yùn)動(dòng)初始狀態(tài)決定，即預(yù)報(bào)偏差矢量隨時(shí)間的演化過程是確定性的。只要獲取多組不同時(shí)刻的相對運(yùn)動(dòng)狀態(tài)量或其線性組合，即可對未知量相對運(yùn)動(dòng)初始狀態(tài)進(jìn)行求解，進(jìn)而得到任意時(shí)刻t的相對運(yùn)動(dòng)狀態(tài)

1.3 角度觀測模型

預(yù)報(bào)偏差矢量在空間中的演化過程對地面觀測站而言表現(xiàn)為預(yù)報(bào)目標(biāo)與真實(shí)目標(biāo)在光學(xué)測角平面上投影點(diǎn)的相對運(yùn)動(dòng)。為了量化這種映射關(guān)系，構(gòu)建如圖2 所示的角度觀測模型。其中，圖2 所示平面為光學(xué)測角平面，O'與P'分別為空間幾何模型中預(yù)報(bào)目標(biāo)O與真實(shí)目標(biāo)P在光學(xué)測角平面上的投影點(diǎn)，v是預(yù)報(bào)目標(biāo)O的空間飛行速度方向在光學(xué)測角平面上的投影。定義α軸方向?yàn)榕cv平行的沿跡方向，β軸方向?yàn)榇怪庇趘的垂跡方向，若將角度偏差觀測量按所定義的沿跡方向和垂跡方向進(jìn)行分解和量化，則P'在以O(shè)'為原點(diǎn)的O'αβ直角坐標(biāo)平面內(nèi)的坐標(biāo)(α，β)即為角度觀測模型的角度偏差觀測量。

聯(lián)合所構(gòu)建的空間幾何模型與角度觀測模型，可得到角度偏差觀測量(α，β)與Oxyz坐標(biāo)系下相對位置狀態(tài)量(x，y，z)的幾何映射關(guān)系為

式中，c1=cosγ1，s1=sinγ1，c2=cosγ2，s2=sinγ2，ρ為距離預(yù)報(bào)值，對應(yīng)于預(yù)報(bào)矢量的模長。同時(shí)可得到距離預(yù)報(bào)偏差Δρ與相對位置狀態(tài)量(x，y，z)的幾何映射關(guān)系為

式（4）和式（5）所描述的是某一時(shí)刻的靜態(tài)幾何映射關(guān)系，而其中c1，s1，c2，s2，ρ均隱含為時(shí)間t的函數(shù)，與角度偏差的觀測時(shí)刻相對應(yīng)，即幾何映射關(guān)系隨時(shí)間動(dòng)態(tài)更新。聯(lián)合式（4）的幾何映射關(guān)系與狀態(tài)演化模型可知，只要獲取多組不同時(shí)刻的角度偏差觀測量，即可映射得到多組相對運(yùn)動(dòng)初始狀態(tài)的線性組合，進(jìn)而實(shí)現(xiàn)對相對運(yùn)動(dòng)初始狀態(tài)的求解。進(jìn)一步地，只要得到相對運(yùn)動(dòng)初始狀態(tài)的近似解或最優(yōu)估計(jì)值，即可通過式（2）得到任意時(shí)刻t的相對位置狀態(tài)量(xt，yt，zt)，并通過式（5）的幾何映射關(guān)系最終得到各時(shí)刻所對應(yīng)的距離預(yù)報(bào)偏差值。

2 算法設(shè)計(jì)

距離預(yù)報(bào)偏差實(shí)時(shí)修正算法設(shè)計(jì)的關(guān)鍵是通過角度偏差觀測量的累積使距離預(yù)報(bào)偏差估計(jì)值迅速向真值收斂。從狀態(tài)演化模型中C-W 方程解的形式可知，其描述的是由相對運(yùn)動(dòng)初始狀態(tài)所定義的線性動(dòng)態(tài)系統(tǒng)，可使用經(jīng)典的卡爾曼濾波算法對其進(jìn)行狀態(tài)估計(jì)。定義X=為卡爾曼濾波器的狀態(tài)向量，Z=(ρα ρβ)T為觀測向量，其無控狀態(tài)方程和觀測方程分別為

式中，wk-1為過程噪聲向量，vk為觀測噪聲向量，wk-1和vk均為白噪聲向量且相互獨(dú)立。式（6）的狀態(tài)方程中Xk和Xk-1均對應(yīng)于相對運(yùn)動(dòng)初始狀態(tài)，因此狀態(tài)轉(zhuǎn)移矩陣Ak取單位矩陣E。定義與式（4）對應(yīng)的幾何映射矩陣H為

與式（2）和式（3）對應(yīng)的相對運(yùn)動(dòng)狀態(tài)轉(zhuǎn)移矩陣M為

則式（7）的觀測方程中觀測矩陣Ck可表示為

即通過幾何映射矩陣H和相對運(yùn)動(dòng)狀態(tài)轉(zhuǎn)移矩陣M，角度偏差觀測量與相對運(yùn)動(dòng)初始狀態(tài)直接建立了映射關(guān)系，其中矩陣H隱含有時(shí)間變量t，且H和M中的時(shí)刻t與角度偏差量的觀測時(shí)刻t=tk相對應(yīng)。

卡爾曼濾波對狀態(tài)量的迭代求解包括預(yù)測和更新兩個(gè)步驟。其預(yù)測過程如式（11），包括狀態(tài)向量的先驗(yàn)估計(jì)與協(xié)方差矩陣計(jì)算。算法設(shè)定狀態(tài)向量估計(jì)值的初值為零向量，同時(shí)設(shè)定狀態(tài)協(xié)方差矩陣Pk的初值P0=diag((Δx)2(Δy)2(Δz)2，其中相對位置標(biāo)準(zhǔn)差Δx，Δy，Δz均取100 m，相對速度標(biāo)準(zhǔn)差均取0.1 m/s，過程噪聲協(xié)方差矩陣Qk取零矩陣。

卡爾曼濾波的更新過程如式（12），包括卡爾曼增益Kk的計(jì)算以及狀態(tài)向量估計(jì)值與協(xié)方差矩陣Pk的更新過程。算法設(shè)定觀測噪聲協(xié)方差矩陣Rk=diag((ρΔα)2(ρΔβ)2)，其中Δα，Δβ均與角度觀測量標(biāo)準(zhǔn)差σA相等，對應(yīng)于望遠(yuǎn)鏡系統(tǒng)的角度測量精度指標(biāo)，設(shè)置為5″以內(nèi)，ρ對應(yīng)于t=tk時(shí)刻的距離預(yù)報(bào)值。

式（11）和式（12）的卡爾曼濾波算法作為最優(yōu)狀態(tài)求解器實(shí)時(shí)迭代產(chǎn)生相對運(yùn)動(dòng)初始狀態(tài)估計(jì)值，在每一次迭代計(jì)算的同時(shí)，將產(chǎn)生的代入式（2）可得到t=tk時(shí)刻的相對位置狀態(tài)量估計(jì)值而后更新式（5）的幾何映射關(guān)系并將代入，即得到了本次迭代對應(yīng)的距離預(yù)報(bào)偏差估計(jì)值。

在算法的實(shí)際應(yīng)用中，卡爾曼濾波器的角度偏差觀測量輸入由空間碎片真實(shí)位置觀測值和預(yù)報(bào)位置之間的角度偏差依所建立的角度觀測模型實(shí)時(shí)轉(zhuǎn)換得到，其中真實(shí)位置通過方位俯仰或赤經(jīng)赤緯形式由觀測系統(tǒng)實(shí)時(shí)提供，與望遠(yuǎn)鏡的指向方向無關(guān)，而各時(shí)刻的預(yù)報(bào)位置在觀測開始之前即已計(jì)算得到，且作為已知量在觀測過程中不再改變。即算法運(yùn)行的距離預(yù)報(bào)偏差修正過程與目標(biāo)跟蹤過程中的望遠(yuǎn)鏡指向偏差修正動(dòng)作相互獨(dú)立，互不影響。

在觀測弧段的起始，系統(tǒng)開始接收光學(xué)測角數(shù)據(jù)并實(shí)時(shí)轉(zhuǎn)換為卡爾曼濾波器的角度偏差觀測量輸入，算法迭代產(chǎn)生的距離預(yù)報(bào)偏差估計(jì)值被實(shí)時(shí)傳輸至測距系統(tǒng)。測距系統(tǒng)依據(jù)接收到的距離預(yù)報(bào)偏差修正量實(shí)時(shí)對TLE 距離預(yù)報(bào)進(jìn)行修正并搜索回波，直至弧段結(jié)束，由此即實(shí)現(xiàn)了距離預(yù)報(bào)偏差實(shí)時(shí)修正技術(shù)輔助的空間碎片測距過程。

3 數(shù)值仿真與結(jié)果分析

為了驗(yàn)證距離預(yù)報(bào)偏差實(shí)時(shí)修正算法的有效性，選取軌道高度在1 500 km 以下的9 顆ILRS 近地聯(lián)測激光衛(wèi)星作為數(shù)值仿真對象，基本涵蓋了近地目標(biāo)不同的軌道高度和軌道傾角。所選目標(biāo)同時(shí)具有TLE 兩行根數(shù)預(yù)報(bào)和ILRS 提供的高精度CPF 預(yù)報(bào)，其中無修正TLE 預(yù)報(bào)軌道與CPF 預(yù)報(bào)軌道分別對應(yīng)算法模型中的預(yù)報(bào)目標(biāo)軌道和真實(shí)目標(biāo)軌道。本文研究的是距離預(yù)報(bào)偏差改進(jìn)問題，將CPF 預(yù)報(bào)軌道作為參考軌道，其軌道精度相比TLE 預(yù)報(bào)能夠滿足算法驗(yàn)證的軌道精度需求［2］。

為了模擬真實(shí)測距場景，在數(shù)值仿真中將CPF 角度預(yù)報(bào)作為卡爾曼濾波器的角度觀測量輸入，同時(shí)將CPF 距離預(yù)報(bào)作為真實(shí)距離觀測值參考。定義CPF 距離預(yù)報(bào)與無修正TLE 距離預(yù)報(bào)之間的距離偏差為參考距離偏差ΔρREF，則算法實(shí)時(shí)產(chǎn)生的距離偏差估計(jì)值向參考距離偏差ΔρREF的收斂速度以及與其接近程度將作為數(shù)值仿真結(jié)果的評價(jià)標(biāo)準(zhǔn)。

獲取9 顆近地目標(biāo)于2021-12-01～2021-12-20 日間在長春站過境的所有弧段，在篩除最大俯仰角小于10°以及初始距離偏差小于100 m 的弧段后共計(jì)633 圈，各弧段取10°作為觀測的起止俯仰角并做預(yù)報(bào)。從實(shí)用化的角度，將TLE 距離預(yù)報(bào)值與真實(shí)距離觀測值之間的偏差修正至小于100 m 即可滿足空間碎片觀測時(shí)段常規(guī)探測器的距離門寬范圍，從而有效增加回波捕獲概率。因此，數(shù)值仿真以從弧段起始至|ΔρREF|＜100 所用時(shí)長作為修正時(shí)間，將修正時(shí)間與修正時(shí)間的全弧段時(shí)長占比作為評價(jià)指標(biāo)。

仿真并統(tǒng)計(jì)各目標(biāo)20 天內(nèi)所有弧段的平均弧段時(shí)長、平均初始距離偏差、平均最大距離偏差、平均修正時(shí)間與平均修正時(shí)間弧段時(shí)長占比，各目標(biāo)仿真統(tǒng)計(jì)結(jié)果如表1 所示，其中弧段時(shí)長和修正時(shí)間單位為分鐘（min）。為了仿真對比角度測量精度對算法修正精度和收斂速度的影響，并驗(yàn)證算法對于角度測量精度指標(biāo)的寬容度，將卡爾曼濾波器中角度觀測量標(biāo)準(zhǔn)差σA分別設(shè)置為2″和5″。

從仿真統(tǒng)計(jì)結(jié)果可以看出，所選近地合作目標(biāo)的平均弧段時(shí)長在16 min 以內(nèi)，平均初始距離偏差量均小于500 m，且弧段最大距離偏差略大于初始距離偏差，兩者相比差別不大。角度觀測量標(biāo)準(zhǔn)差σA為2″時(shí)，算法能夠在大部分弧段起始的1 min 以內(nèi)將TLE 預(yù)報(bào)的距離偏差量修正至小于100 m，各目標(biāo)平均修正時(shí)間弧段時(shí)長占比最大值小于15%，即大部分弧段能夠在弧段起始的前15%完成距離偏差修正過程。對比σA為5″的仿真結(jié)果，平均修正時(shí)間較σA為2″時(shí)略有增加，表明較高的角度測量精度有利于加快距離偏差的修正速度，各目標(biāo)平均修正時(shí)間弧段時(shí)長占比最大值小于20%，仍在可接受范圍內(nèi)，即算法對于角度測量精度指標(biāo)具有較好的寬容度，能夠適用于具備不同角度測量精度的望遠(yuǎn)鏡系統(tǒng)。

以所選目標(biāo)中平均初始距離偏差最大的GRACE-FO-1 為例，圖3 顯示了GRACE-FO-1 于協(xié)調(diào)世界時(shí)（Universal Time Coordinated，UTC）2021 年12 月18 日01：29：17～01：36：21 在長春站過境弧段的完整距離預(yù)報(bào)偏差修正過程（σA=2″）。

圖3 中弧段時(shí)長7.1 min，初始距離偏差接近1 000 m。在弧段的起始，距離偏差估計(jì)值的初值為0，隨著光學(xué)角度偏差觀測量的累積，修正算法使得距離偏差估計(jì)值迅速向參考距離偏差ΔρREF收斂；在弧段開始1 min 后，ΔρREF與的偏差量減小至100 m 以內(nèi)，也即TLE 距離預(yù)報(bào)值與真實(shí)距離觀測值之間的偏差量被修正至100 m 以內(nèi)。此時(shí)已滿足了常規(guī)探測器的距離門寬搜索范圍要求，激光測距系統(tǒng)可穩(wěn)定地獲取有效回波，隨后距離偏差量進(jìn)一步減小，一直穩(wěn)定在ΔρREF附近，直至弧段結(jié)束。

本文算法主要在于通過減小距離預(yù)報(bào)偏差來縮減距離門回波搜索范圍，提高有效回波信號(hào)落入距離門寬內(nèi)的概率，進(jìn)而提升距離門控回波搜索捕獲效率和整體測距成功概率。其對于測距成功概率的提升程度則主要與不同觀測時(shí)段不同天空背景噪聲強(qiáng)度下所允許的探測器最大距離門寬，以及距離預(yù)報(bào)偏差修正前后的距離偏差量有關(guān)?？臻g碎片測距能力的提升是多種技術(shù)手段綜合運(yùn)用的結(jié)果，在實(shí)際應(yīng)用中，將本文提出的距離預(yù)報(bào)偏差實(shí)時(shí)修正技術(shù)與其它測距技術(shù)手段相結(jié)合，能夠發(fā)揮出更大的作用。

4 結(jié)論

非合作目標(biāo)空間碎片的距離預(yù)報(bào)偏差修正問題一直是領(lǐng)域內(nèi)的難題，將視位置角度偏差作為觀測量的距離預(yù)報(bào)偏差實(shí)時(shí)修正方法為其提供了新的思路和解決方案，仿真統(tǒng)計(jì)結(jié)果表明方法具有很高的實(shí)用價(jià)值，而能夠持續(xù)穩(wěn)定地獲取角度偏差觀測量是方法成功的關(guān)鍵。當(dāng)前空間碎片白天測距仍是領(lǐng)域內(nèi)的難題，受限于白天強(qiáng)天空背景噪聲下的探測器最大距離門寬限制，其對于距離預(yù)報(bào)精度的要求更加嚴(yán)苛，若能夠很好地解決白天空間目標(biāo)的可視問題，方法將在空間碎片白天測距中也發(fā)揮重要的作用，有利于空間碎片白天測距能力的提升。算法模型中指向偏差、蒙氣差、角度采樣延時(shí)、大氣延遲誤差等觀測量和狀態(tài)量噪聲對于算法穩(wěn)定性的影響將通過進(jìn)一步的數(shù)值仿真和實(shí)測實(shí)驗(yàn)進(jìn)行研究和驗(yàn)證。