廖祥, 鄭靖, 許永生, 謝成清, 丁強(qiáng)強(qiáng)

(1. 深圳市魔方衛(wèi)星科技有限公司, 深圳 518129; 2. 深圳航天東方紅衛(wèi)星有限公司, 深圳 518057)

2015 年以來(lái),多個(gè)國(guó)外商業(yè)航天初創(chuàng)公司相繼提出并部署各自的輕小型合成孔徑雷達(dá)(mini synthetic aperture radar,miniSAR)衛(wèi) 星[1]。 在 此背景下,深圳航天東方紅衛(wèi)星有限公司瞄準(zhǔn)市場(chǎng)前景廣闊的高分辨率微波成像和形變測(cè)量應(yīng)用領(lǐng)域的miniSAR 衛(wèi)星項(xiàng)目,并立項(xiàng)研發(fā),計(jì)劃2021年底具備發(fā)射條件。 為實(shí)現(xiàn)高精度長(zhǎng)時(shí)間序列的形變監(jiān)測(cè),miniSAR 衛(wèi)星的重軌干涉基線設(shè)計(jì)優(yōu)于500 m,要求衛(wèi)星在軌任務(wù)期間的三維空間回歸軌跡被控制在直徑為500 m 的管道內(nèi)。

管道指以衛(wèi)星運(yùn)行軌道的參考軌跡為中心、干涉基線極大值為直徑的三維空間[2]。 以國(guó)外早期干涉SAR 衛(wèi)星項(xiàng)目為例,ERS-1/ERS-2 管道設(shè)計(jì)直徑為600 m[3],TerraSAR-X 為500 m[4],ALOS/ALOS-2 為1000 m[5],Sentinel-1 則 達(dá) 到 了150 m[6-7]。 由于SAR 衛(wèi)星一般運(yùn)行在500 ～800 km高度范圍內(nèi)的低地球軌道空間,大氣阻力與太陽(yáng)光壓等非保守?cái)z動(dòng)力會(huì)對(duì)衛(wèi)星軌道根數(shù)造成長(zhǎng)期影響,如此高精度的管道控制要求,需設(shè)計(jì)高精度參考軌跡為基準(zhǔn),進(jìn)一步實(shí)施精密軌道控制[8]。而參考軌跡由參考軌道預(yù)報(bào)一個(gè)回歸周期生成,因此,干涉SAR 衛(wèi)星高精度軌道控制中一個(gè)十分重要的術(shù)語(yǔ)“參考軌道”被提出。

參考軌道指在僅考慮中心天體保守力情況下,滿足預(yù)報(bào)一個(gè)回歸周期后衛(wèi)星地固系下的三維空間位置完成重合的一類軌道。 參考軌道僅為一組包含歷元的軌道根數(shù),其最終類型為瞬時(shí)軌道根數(shù)(又稱密切軌道根數(shù)),坐標(biāo)系通常為慣性系瞬時(shí)真赤道地心坐標(biāo)系(true of date,TOD)[9]。對(duì)于需要生成長(zhǎng)期一致的時(shí)間序列干涉數(shù)據(jù)的SAR 衛(wèi)星,參考軌道一旦完成設(shè)計(jì),在整個(gè)任務(wù)周期內(nèi)就不能更改,這同時(shí)對(duì)重軌單星和雙星近距編隊(duì)干涉SAR 系統(tǒng)構(gòu)成約束。

由于干涉SAR 衛(wèi)星技術(shù)起步較晚,且工程任務(wù)較少,國(guó)內(nèi)外關(guān)于參考軌道設(shè)計(jì)的文獻(xiàn)較欠缺。推導(dǎo)嚴(yán)密、可復(fù)現(xiàn)性好的精密參考軌道設(shè)計(jì)文獻(xiàn)更是稀缺。 文獻(xiàn)[10]針對(duì)TerraSAR-X 任務(wù),實(shí)現(xiàn)了一種基于虛擬速度增量與序列二次規(guī)劃的參考軌道優(yōu)化設(shè)計(jì)方法,但由于引入了2 次虛擬軌控,生成的參考軌跡不連續(xù);文獻(xiàn)[11]提出基于迭代修正方法的嚴(yán)格回歸軌道設(shè)計(jì),但文中的推導(dǎo)存在部分錯(cuò)誤,實(shí)現(xiàn)算法存在缺陷,且回歸精度較差;文獻(xiàn)[12]提出一種基于高階地球重力場(chǎng)模型的嚴(yán)格回歸軌道單目標(biāo)優(yōu)化設(shè)計(jì)方法,但回歸精度較差,仍在米級(jí);文獻(xiàn)[13]在修正文獻(xiàn)[10]方法的同時(shí),解決了利用傾角偏置修正沿跡向長(zhǎng)期漂移的問(wèn)題,但文中采用專有工具優(yōu)化設(shè)計(jì)參考軌道,無(wú)法獲取具體實(shí)現(xiàn)算法。

本文基于嵌套式迭代修正思想,通過(guò)合理選擇參考系、初始預(yù)估參數(shù)、平根與軌道外推模型,并合理設(shè)置迭代修正反饋值及迭代收斂判據(jù),建立了一套推導(dǎo)嚴(yán)密、可復(fù)現(xiàn)性好的miniSAR 衛(wèi)星精密參考軌道設(shè)計(jì)算法,由其設(shè)計(jì)的參考軌道生成的參考軌跡在地固系下的閉環(huán)精度達(dá)毫米級(jí)。

1 參考軌道特性分析

miniSAR 衛(wèi)星的參考軌道必須同時(shí)具備2 個(gè)基本特性:①軌道必須是回歸軌道,即運(yùn)行在該軌道上的衛(wèi)星對(duì)星下點(diǎn)具有高精度重訪特性;②軌道必須是凍結(jié)軌道,即衛(wèi)星星下點(diǎn)重訪時(shí)軌道高度的一致性。 根據(jù)軌道動(dòng)力學(xué)原理,回歸軌道特性與軌道半長(zhǎng)軸a和傾角i相對(duì)應(yīng),凍結(jié)軌道特性與偏心率e和近地點(diǎn)幅角ω相對(duì)應(yīng)[10]。

此外,考慮到衛(wèi)星能源的穩(wěn)定性及升交點(diǎn)赤經(jīng)的進(jìn)動(dòng)速率,miniSAR 衛(wèi)星的參考軌道還設(shè)計(jì)具備太陽(yáng)同步軌道特性。 這一特性進(jìn)一步約束了參考軌道半長(zhǎng)軸a和傾角i之間的關(guān)系[14]。

1.1 太陽(yáng)同步特性

太陽(yáng)同步軌道平面法線與太陽(yáng)矢量在赤道平面上的投影保持著固定的夾角關(guān)系,這使得運(yùn)行其上的衛(wèi)星可獲得穩(wěn)定的光照條件,且在衛(wèi)星飛越相同地理緯度的目標(biāo)區(qū)域處于同一地方時(shí)。 太陽(yáng)同步軌道的設(shè)計(jì)主要考慮降交點(diǎn)平地方時(shí)(mean local time, MLT)與平均軌道高度(mean orbit altitude, MOA)2 個(gè)參數(shù)。

文獻(xiàn)[15]指出,太陽(yáng)同步軌道傾角平根受太陽(yáng)引力影響的長(zhǎng)期攝動(dòng)年變化率di/dt(°/年)與降交點(diǎn)MLT 成正弦關(guān)系,降交點(diǎn)MLT 為3:00、9:00、15:00、21:00 時(shí)di/dt最 大,降 交 點(diǎn)MLT 為0:00、6:00、12:00、18:00 時(shí)di/dt最小并接近于0,如圖1 所示。 因此,在設(shè)計(jì)太陽(yáng)同步軌道參數(shù)時(shí),干涉SAR 衛(wèi)星適合選擇0:00、6:00、12:00、18:00 附近的降交點(diǎn)MLT,以減小傾角平根的長(zhǎng)期漂移,進(jìn)而避免不必要的傾角調(diào)整。

圖1 受太陽(yáng)引力影響的傾角長(zhǎng)期項(xiàng)攝動(dòng)[15]Fig.1 Secular perturbation of inclination due to gravity of the Sun[15]

MOA 的設(shè)計(jì)則需考慮SAR 衛(wèi)星載荷在一個(gè)回歸周期T內(nèi)可實(shí)現(xiàn)對(duì)全球的無(wú)縫覆蓋。 這里的覆蓋能力并非實(shí)際可照射區(qū)域,而是衛(wèi)星在姿態(tài)最大側(cè)擺角下載荷所能覆蓋的區(qū)域來(lái)計(jì)算。 由于SAR 載荷一般具有照射盲區(qū),因而MOA 的設(shè)計(jì)需要對(duì)此因素加以建模分析。

確定好降交點(diǎn)MLT 與MOA 后,可根據(jù)基于J2低階引力攝動(dòng)下的太陽(yáng)同步軌道設(shè)計(jì)方程進(jìn)一步確定軌道傾角:

式中:μ為地球引力常數(shù);R⊕為地球平均半徑;為地球公轉(zhuǎn)升交點(diǎn)赤經(jīng)進(jìn)動(dòng)速率;J2為地球非球形引力模型二階帶諧項(xiàng)系數(shù)。

1.2 回歸特性

理想的參考軌跡在地固系三維空間下為首尾完全閉合的曲線,因而參考軌道必須具有回歸特性。 給定回歸周期T與對(duì)應(yīng)的回歸圈數(shù)N,則回歸軌道約束方程為

式中:ω⊕為地球自轉(zhuǎn)速率;分別為升交點(diǎn)赤經(jīng)、近地點(diǎn)幅角、平近點(diǎn)角的進(jìn)動(dòng)速率,三者在J2低階引力攝動(dòng)下的解析式[16]為

式中:傾角i以式(1)為約束參與式(3)對(duì)半長(zhǎng)軸a的求解;Ω為升交點(diǎn)赤經(jīng);ω為近地點(diǎn)幅角;M為平近點(diǎn)角。

1.3 凍結(jié)特性

理想?yún)⒖架壽E的星下點(diǎn)重訪時(shí)軌道高度保持不變,因而參考軌道必須具有凍結(jié)特性。 凍結(jié)軌道有2 種類型:①臨界傾角的凍結(jié)軌道;②小偏心率近圓凍結(jié)軌道。 干涉SAR 衛(wèi)星的參考軌道屬于第2 種情況。 凍結(jié)軌道設(shè)計(jì)方程為

根據(jù)小偏心率近圓凍結(jié)軌道的要求,考慮J2及J3低階引力攝動(dòng)可解算得到

式中:J3為地球非球形引力模型三階帶諧項(xiàng)系數(shù)。

1.4 參考軌道預(yù)估值

綜上對(duì)參考軌道特性的分析,可預(yù)估得到5 個(gè)基于低階引力勢(shì)場(chǎng)的參考軌道根數(shù):

需要指明的是,此處獲得的參考軌道初值為T(mén)OD 坐標(biāo)系下的平根,在使用之前須將其轉(zhuǎn)換為瞬根方可作為初始輸入?yún)⑴c后續(xù)迭代修正過(guò)程:

式中:a0、i0、e0、ω0、M0分別為半長(zhǎng)軸、傾角、偏心率、近地點(diǎn)幅角、平近點(diǎn)角的瞬根。

2 參考軌道設(shè)計(jì)算法

參考軌道的設(shè)計(jì)包含兩部分工作:①針對(duì)(a,i)組合的太陽(yáng)同步回歸軌道的優(yōu)化設(shè)計(jì);②針對(duì)(e,ω)組合的凍結(jié)軌道的優(yōu)化設(shè)計(jì)。 但由于(a,i)迭代修正的輸入為預(yù)估的參考軌道根數(shù),2種組合優(yōu)化設(shè)計(jì)方法需嵌套使用。 當(dāng)(e,ω)迭代收斂之后,只需再重復(fù)一輪(a,i)迭代修正即可實(shí)現(xiàn)高精度回歸。 具體算法流程如圖2 所示。

圖2 精密參考軌道設(shè)計(jì)算法流程Fig.2 Algorithm flow chart of precise reference orbit design

2.1 迭代修正具體實(shí)現(xiàn)

從算法流程可見(jiàn),參考軌道設(shè)計(jì)過(guò)程主要分為3 個(gè)步驟:①(a,i)組合迭代粗修正;②(e,ω)組合迭代修正;③(a,i)組合迭代精修正。

2.1.1 (a,i)組合迭代粗修正

(a,i)組合迭代粗修正以參考軌道初值a0為輸入,外推1 個(gè)回歸周期T,利用首尾軌道數(shù)據(jù)生成平根修正量(Δˉa,Δˉi),并進(jìn)一步通過(guò)平瞬轉(zhuǎn)換獲得新的(a,i)值,其余軌道根數(shù)保持不變。 其中,(ˉ)為實(shí)際獲得的平根修正量,(a,i)則為最終修正得到的瞬根值。 ()根據(jù)梯度矩陣獲得

式中:Δλ和Δφ分別為首尾軌道數(shù)據(jù)星下點(diǎn)經(jīng)緯度之差;4 個(gè)偏導(dǎo)數(shù)求解式為

其中:中間量?u/?a計(jì)算如下:

式(8) ～式(10)為最終算法,其詳細(xì)推導(dǎo)過(guò)程可參考文獻(xiàn)[10-12],但注意式(9)對(duì)文獻(xiàn)[11]中錯(cuò)誤部分的更正。

2.1.2 (e,ω)組合迭代修正

圖3 偏心率平根矢量修正示意圖Fig.3 Diagram of mean eccentricity vector correction

2.1.3 (a,i)組合迭代精修正

(a,i)組合迭代精修正的實(shí)施方法與2.2 節(jié)

步驟1相同,但本節(jié)迭代修正的收斂判據(jù)與2.2 節(jié)步驟1 不同。

2.2 收斂判據(jù)

(a,i)與(e,ω)組合每輪修正需要多次迭代實(shí)現(xiàn)收斂,參考軌道設(shè)計(jì)算法的3 輪組合修正都需要給定明確的收斂判據(jù),用以判斷每次迭代后的結(jié)果是否滿足收斂條件,各輪修正的收斂判據(jù)如下。

步驟1判定3。 以前后2 次迭代修正后的參考軌道根數(shù),預(yù)報(bào)1 個(gè)回歸周期T之后的三維空間軌跡回歸精度優(yōu)于精密定軌精度為收斂準(zhǔn)則。 參考GNSS 精密定軌精度為厘米級(jí),故收斂準(zhǔn)則設(shè)置為

式中:n為迭代次數(shù),n≥1;為參考軌跡閉環(huán)回歸精度,計(jì)算式為

式中:xf、yf、zf分別為地固系三維坐標(biāo)值,時(shí)間下標(biāo)t0、t0+T表征預(yù)報(bào)1 個(gè)回歸周期T所得軌道數(shù)據(jù)的首、尾坐標(biāo)值。

步驟2判定2。 以前后2 次迭代修正后的偏心率平根矢量偏差引起的參考軌跡偏差ˉaΔˉe優(yōu)于厘米級(jí)為收斂準(zhǔn)則。 由于ˉa為10-6m量級(jí),故偏心率平根矢量前后差異需優(yōu)于1.0 ×10-8,即

步驟3判定1。 以前后2 次迭代修正后的參考軌道根數(shù)預(yù)報(bào)1 個(gè)回歸周期T之后的三維空間軌跡回歸精度無(wú)明顯提高為收斂準(zhǔn)則。 考慮二階展開(kāi),回歸精度差異需優(yōu)于o(),即

式中:abs()表示取絕對(duì)值;n為迭代次數(shù),n≥2。

2.3 平瞬轉(zhuǎn)換

參考軌道設(shè)計(jì)算法的迭代實(shí)際修正量為平根,但參考軌道的形式為瞬根,因此在修正平根后需執(zhí)行平瞬轉(zhuǎn)換獲得瞬根修正值。

常用的平瞬轉(zhuǎn)換模型有多種[16],本文選用考慮J2～J6及攝動(dòng)的Eckstein-Hechler 模型[17],該模型專門(mén)針對(duì)小偏心率凍結(jié)軌道進(jìn)行優(yōu)化,尤其適合參考軌道設(shè)計(jì)過(guò)程。

Eckstein-Hechler 模型在基于SCILAB 科學(xué)計(jì)算語(yǔ)言的Celestlab[18]與基于Java 語(yǔ)言的Orekit[19]這2 套成熟開(kāi)源航天動(dòng)力學(xué)庫(kù)中皆有實(shí)現(xiàn),可作為實(shí)現(xiàn)參考。

2.4 外推模型

(a,i)與(e,ω)組合迭代修正的過(guò)程皆需要進(jìn)行軌道外推,以獲得向后一段時(shí)間的軌道數(shù)據(jù)。由于參考軌道只考慮中心天體保守力,因而外推模型的動(dòng)力學(xué)只考慮高階地球非球形引力,本文為100 ×100 階的EGM2008 引力模型[20]。

動(dòng)力學(xué)模型的單一決定了外推模型的精度取決于積分器的選擇及對(duì)積分步長(zhǎng)的控制。 考慮到積分效率與精度,選擇9 階誤差控制的8 階Runge-Kutta-Verner 積分器RKV8(9)[21]。 修正過(guò)程中,誤差容限設(shè)為優(yōu)于1.0 ×10-13,(a,i)、(e,ω)組合迭代修正最大步長(zhǎng)分別控制在10 s、5 s 以內(nèi)。

3 數(shù)值實(shí)驗(yàn)

3.1 輸入數(shù)據(jù)與實(shí)驗(yàn)條件

3.1.1 軌道特性輸入

miniSAR 衛(wèi)星軌道設(shè)計(jì)為具有10 天/151 軌(回歸周期為10 天,10 天內(nèi)衛(wèi)星運(yùn)行151 軌)回歸、凍結(jié)特性的太陽(yáng)同步軌道,標(biāo)稱降交點(diǎn)MLT為6:00(晨昏軌道),衛(wèi)星質(zhì)量為180 kg。

3.1.2 定位數(shù)據(jù)輸入

衛(wèi)星成功入軌后經(jīng)運(yùn)控系統(tǒng)提供一組包含升交點(diǎn)的定位數(shù)據(jù),如表1 所示(表中為仿真的地固系三維位置與速度值)。 對(duì)其插值得到升交點(diǎn)處地固系狀態(tài)矢量,對(duì)應(yīng)歷元t0為2021-01-01,00:26:13.511 668 UTCG。 再利用升交點(diǎn)處的地固狀態(tài)矢量計(jì)算得到升交點(diǎn)赤經(jīng)平根初值為=10.881 803 934°。

表1 定位數(shù)據(jù)Table 1 Precision positioning data

3.1.3 參考軌道初值

根據(jù)軌道特性預(yù)估值與定位數(shù)據(jù),計(jì)算得到的參考軌道在t0時(shí)刻的初值為

平瞬轉(zhuǎn)換后參與迭代修正的參考軌道初值為

3.1.4 實(shí)驗(yàn)條件

數(shù)值實(shí)驗(yàn)基于MATLAB 語(yǔ)言環(huán)境,軌道外推與數(shù)據(jù)報(bào)表使用STK Engine 的MATLAB 編程接口實(shí)現(xiàn),平瞬轉(zhuǎn)換算法移植自Celestlab。 仿真程序運(yùn)行于搭載Windows 10 專業(yè)版操作系統(tǒng)的聯(lián)想P52 移動(dòng)工作站,機(jī)帶內(nèi)存32 GB,處理器型號(hào)為Intel Xeon E-2176M。

3.2 實(shí)驗(yàn)結(jié)果

表2 給出了參考軌道初值及數(shù)值實(shí)驗(yàn)過(guò)程中3 輪迭代修正所得參考軌道根數(shù),其中a為最終獲得的參考軌道根數(shù)。 從表2 中數(shù)據(jù)可見(jiàn),由a0到a1及a2到a的過(guò)程中,只有半長(zhǎng)軸與傾角是變化量;由a1到a2的過(guò)程中,只有偏心率、近地點(diǎn)幅角與平近點(diǎn)角是變化量。

表2 迭代修正過(guò)程中所得參考軌道根數(shù)(TOD 瞬根)Table 2 Reference orbit elements obtained during process of iterative correction (TOD osculating elements)

表3 給出了3 輪修正迭代過(guò)程中收斂判據(jù)的統(tǒng)計(jì)數(shù)據(jù),圖4 ～圖7 則對(duì)統(tǒng)計(jì)數(shù)據(jù)進(jìn)行了可視化。 其中,圖4 與圖5 對(duì)應(yīng)第1、2 輪(a,i)組合迭代修正過(guò)程中的收斂情況,橫坐標(biāo)n為迭代次數(shù),圖6 與圖7 則對(duì)應(yīng)(e,ω)組合迭代修正過(guò)程中偏心率平根矢量收斂情況。 統(tǒng)計(jì)數(shù)據(jù)表明:

圖4 第1 輪(a,i) 迭代修正收斂情況Fig.4 Convergence of first (a,i) iterative correction

圖5 第2 輪(a,i) 迭代修正收斂情況Fig.5 Convergence of second (a,i) iterative correction

圖6 偏心率平根矢量收斂情況Fig.6 Convergence of mean eccentricity vector

圖7 (e,ω)迭代修正收斂情況Fig.7 Convergence of (e,ω) iterative correction

表3 迭代修正收斂數(shù)據(jù)Table 3 Convergence data of iterative correction

1) 第1 輪(a,i)組合迭代修正過(guò)程中,僅6 次迭代即實(shí)現(xiàn)收斂,收斂精度達(dá)4.38 ×10-7m。

2) (e,ω)組合迭代修正過(guò)程中,120 天內(nèi)以1 個(gè)回歸周期為步長(zhǎng)的偏心率平根矢量在迭代修正過(guò)程中逐步收縮半徑,僅4 次迭代偏心率平根變化率Δˉe即收斂至2 ×10-10。

3) 第2 輪(a,i)組合迭代修正過(guò)程中,僅4 次迭代即實(shí)現(xiàn)收斂,收斂精度優(yōu)于0.006 m。

綜合本節(jié)所述,實(shí)驗(yàn)結(jié)果的收斂精度滿足收斂判據(jù)指標(biāo),驗(yàn)證了所提miniSAR 衛(wèi)星參考軌道迭代修正設(shè)計(jì)算法的有效性和精密性。

實(shí)驗(yàn)中,單次(a,i)修正平均耗時(shí)56.64 s,單次(e,ω)修正平均耗時(shí)624.45 s,合計(jì)耗時(shí)47.29 min。衛(wèi)星軌控任務(wù)周期一般大于軌道周期,合計(jì)耗時(shí)小于miniSAR 衛(wèi)星約95.36 min 的軌道周期,不影響軌控任務(wù)正常開(kāi)展。 此外,考慮到軟硬件條件,參考軌道設(shè)計(jì)工作不適合在內(nèi)存低于GB 級(jí)的嵌入式星載計(jì)算機(jī)上進(jìn)行。

4 結(jié) 論

本文建立了一套完整的miniSAR 衛(wèi)星精密參考軌道設(shè)計(jì)算法,經(jīng)過(guò)嵌套式(a,i)、(e,ω)組合的3 輪迭代修正獲得了毫米級(jí)三維空間軌跡回歸精度。

1) 引入了參考軌道的具體定義,提出的參考軌道設(shè)計(jì)算法流程清晰、推導(dǎo)嚴(yán)密、收斂判據(jù)明確、可復(fù)現(xiàn)性好。

2) 提升了參考軌道設(shè)計(jì)速度與精度,數(shù)值實(shí)驗(yàn)在48 min 內(nèi)獲得了地固系三維空間閉環(huán)回歸精度約為0.005 m 優(yōu)于0.01 m 的精密參考軌道,滿足實(shí)際工程應(yīng)用需求。