林厚源

(1 中國科學(xué)院紫金山天文臺 南京 210023)

(2 中國科學(xué)院空間目標與碎片觀測重點實驗室 南京 210023)

1 引言

北斗G2衛(wèi)星(BeiDou-G2/Compass G2)是北斗二代導(dǎo)航系統(tǒng)的第2顆衛(wèi)星,也是北斗二代中第1顆位于同步軌道(Geosynchronous Orbit,GEO)的G系列衛(wèi)星.它與其他北斗一、二代導(dǎo)航衛(wèi)星一樣是基于東方紅三號衛(wèi)星平臺,其兩側(cè)搭載了太陽能帆板.它于2009年由長征三號丙火箭(CZ-3C)從西昌衛(wèi)星發(fā)射中心發(fā)射上天,但1 yr后便失效,隨后在同步軌道帶(GEO-Belt)[1]開始了大范圍漂移.GEO是氣象、通信和地面監(jiān)測的重要空間,一個漂移的衛(wèi)星對所有在GEO上的衛(wèi)星都是巨大的威脅.為此,在2022年1月,北斗G2被實踐21號衛(wèi)星捕獲并拖進了墳?zāi)管壍?實踐21號衛(wèi)星發(fā)射于2021年,主要用于空間碎片減緩技術(shù)試驗驗證.這次行動展示了我國繼Northrop Grumman公司的Mission Extension Vehicle (MEV)系列衛(wèi)星1之后同樣具有近距離操作、對接和機動的能力,從而能夠?qū)崿F(xiàn)空間碎片的有效清除.但與MEV和衛(wèi)星Intelsat的對接不同的是,北斗G2是顆廢棄的衛(wèi)星.為了能夠安全地實現(xiàn)捕獲,必須要提前知道該捕獲對象的旋轉(zhuǎn)狀態(tài),從而采取正確的策略和操作.當前學(xué)者們還掙扎于如何擬合和反演空間碎片的旋轉(zhuǎn)狀態(tài),尚未有人有效實現(xiàn)空間碎片的旋轉(zhuǎn)狀態(tài)外推計算,因此北斗G2衛(wèi)星的旋轉(zhuǎn)狀態(tài)預(yù)估是個很大的挑戰(zhàn).

在2010-2019年的10 yr里,基于由中國科學(xué)院紫金山天文臺牽頭運行的中國最大的空間碎片光學(xué)觀測網(wǎng),我們采集到了北斗G2衛(wèi)星共計2964個弧段的光度數(shù)據(jù),其中922條可提取出可靠的自轉(zhuǎn)頻率.此外,Graz衛(wèi)星激光測距站額外提供了15條單光子計數(shù)器的光度觀測結(jié)果[2],填補了2015年5月到10月期間由于衛(wèi)星星下點向西漂移超出中國觀測范圍而出現(xiàn)的數(shù)據(jù)空白.我們期望基于這些數(shù)據(jù)建立足夠高精度的旋轉(zhuǎn)運動模型,從而實現(xiàn)北斗G2衛(wèi)星未來旋轉(zhuǎn)狀態(tài)的預(yù)測.

2 北斗G2衛(wèi)星姿軌演化的不同階段

圖1展示了北斗G2衛(wèi)星過去10 yr里自轉(zhuǎn)速度的變化.從圖1(a)可以看出其轉(zhuǎn)速呈振蕩上升的趨勢,而在轉(zhuǎn)速的每月變化趨勢(圖1(b))中可以發(fā)現(xiàn),北斗G2衛(wèi)星轉(zhuǎn)速的振蕩存在明顯的周年效應(yīng),其中12-5月為減速階段,6-11月為加速階段.可推測太陽輻射作用是引起這種變化的主要因素.對此我們以輻射力矩模型[3]為核心建立了動力學(xué)模型并進行數(shù)值擬合和動力學(xué)參數(shù)求解[4].在擬合中我們評估了各種額外的可能引起轉(zhuǎn)速變化的作用因素,見表1,其中只有重力梯度力矩和中心衛(wèi)星主體受到的力矩需要被考慮.由于無法找到滿足整個10 yr演化的統(tǒng)一參數(shù),因此我們采用分段擬合,最終把可擬合的數(shù)據(jù)段命名為S1-S6,其中S6段由于時間跨度太長影響計算效率,被分為S6a和S6b分別擬合(見圖1(c)-(f)).而其他無法被擬合的異常段,結(jié)合軌道根數(shù),可以確定共計6次異常變化,命名為A1-A6 (見圖1(a)).由于錯誤的分段無法使擬合殘差降低到噪聲水平[4],從而確保了分段的可靠性.

表1 數(shù)值模擬中各因素對北斗G2衛(wèi)星轉(zhuǎn)速影響的相對變化量級Table 1 The order of magnitude of the relative deviation caused by each factor for BeiDou-G2 in its variation of rotational speed from the numerical simulation

這些無法由數(shù)值模型擬合的異常段顯示了北斗G2在過去10 yr里過得并不安穩(wěn),這會使得未來旋轉(zhuǎn)狀態(tài)的預(yù)測成為很大的挑戰(zhàn),因此我們必須弄清楚曾經(jīng)發(fā)生了什么.在圖1(a)中將自轉(zhuǎn)頻率數(shù)據(jù)與軌道半長徑和偏心率重疊,可以看到在異常段中,除了A1以外可以分為兩類: A3、A4和A6持續(xù)大約100 d,自轉(zhuǎn)頻率快速增加,同時偏心率也大幅增加,而半長徑幾乎沒有異常變化; 而A2和A5持續(xù)時間大約1個月,期間半長徑發(fā)生過抖動,而偏心率和自轉(zhuǎn)頻率沒有明顯異常.這些異常變化顯然是存在額外且持續(xù)的動力作用,其中燃料泄漏是最可能的因素之一.而這兩不同的變化類型疑似是來自不同的因素,如不同部位的泄漏,或燃料、氣囊等不同類型的泄漏.

圖1 北斗G2衛(wèi)星自轉(zhuǎn)速率過去10 yr的歷史演化和擬合與外推結(jié)果.(a)北斗G2衛(wèi)星自轉(zhuǎn)速度(藍紅,右軸)以及半長徑(青,左一軸)、偏心率(灰,左二軸)變化趨勢.自轉(zhuǎn)數(shù)據(jù)從紫金山天文臺和Graz觀測的光度數(shù)據(jù)中提取,軌道數(shù)據(jù)源于TLE.可擬合的正常轉(zhuǎn)速被分為了S1-S6b共7段用藍色標記,而其中間隔的紅色標記的A1-A6為轉(zhuǎn)速或軌道根數(shù)異常段.據(jù)推斷,在A1處發(fā)生了一次碰撞,在A4段發(fā)生了帆板損壞,而在A6段的起始處發(fā)生了解體事件.(b)北斗G2衛(wèi)星自轉(zhuǎn)速度每月變化趨勢.下面4圖為擬合以及外推結(jié)果,其中(c) S1和S2段; (d) S3,S4和S5段; (e) S6a段; (f)S6b段.各圖底部為帶權(quán)重的擬合殘差.(c)圖中的紫色標記為將S1和S2旋轉(zhuǎn)狀態(tài)的擬合結(jié)果外推到56128 d+12.05 h時刻的結(jié)果.Fig.1 Evolution of the rotation speed of the BeiDou-G2 satellite in the past 10 years and the fitting and propagating results.(a) The BeiDou-G2 satellite’s rotation speed (blue-red,right ordinate),semi-major axis (green,left,first ordinate),and eccentricity (gray,left,second ordinate).The rotation data were extracted from the photometric data that Purple Mountain Observatory and Graz observed,and the orbit data were derived from the two-line element set.The blue dots,which represent the normal rotational speeds that could be fitted,are divided into seven segments,S1-S6b; and the red dots named A1-A6 are in the intervals with abnormal rotational speeds or orbital elements.It was inferred that a collision was included in A1,one solar panel was damaged in A4,and fragmentation occurred at the beginning of A6.(b) Monthly changes of the rotation speed of the BeiDou-G2 satellite.The lower four panels show the results of fitting and propagation in (c) segments S1 and S2; (d)segments S3,S4,and S5; (e) segment S6a; and (f) segment S6b.The bottoms of the four panels are the fitting residuals with weights.The purple in panel (c) of this figure marks the time when the fitting results of the rotation states of S1 and S2 were propagated to 56128 d+12.05 h.

3 自轉(zhuǎn)加速的導(dǎo)火索: 碰撞?

A1段的變化與其他異常段都不相同.北斗G2的偏心率在A1處的變化很突然.雖然缺少當時的轉(zhuǎn)速數(shù)據(jù),但S1和S2的旋轉(zhuǎn)狀態(tài)不一樣(即無法合并在一起擬合).這種變化理論上可能是一次控制機動.但是從結(jié)果來看軌道和姿態(tài)沒有發(fā)生本質(zhì)變化,從時間來看離它最后一次機動控制(官方公布失效)已經(jīng)過去了2 yr,而且前后也沒有類似事件,最關(guān)鍵一點是A1段時候衛(wèi)星正處于非洲上空,超出了中國的可視范圍,沒有理由選在這個時間進行一次主動控制.因此我們深度懷疑這是遭受了一次撞擊.我們通過MJD在56128-56132時間內(nèi)對前后兩組雙行根數(shù)(Two-Line Element,TLE)進行軌道外推,確定了撞擊時刻為56128 d+12.05 h,最近距離為1.5 km,該距離差在高軌目標誤差范圍內(nèi),同時可以得到速度變化Δv=(0.764,-0.628,0.171)m·s-1.由于衛(wèi)星主要部件都為剛體,該碰撞可以看作完全彈性碰撞,則基于同步軌道目標的分布可以估算得到撞擊物的質(zhì)量大概率是一個小于10 kg的空間碎片[4].將S1和S2旋轉(zhuǎn)狀態(tài)的擬合結(jié)果外推到A1位置,則可以進一步確認撞擊細節(jié)[4].撞擊部位可能位于頭部或者尾部(對稱解).

從圖1(a)-(b)中觀察發(fā)現(xiàn),異常事件都發(fā)生在自轉(zhuǎn)加速段(從6月到9月開始),而且自撞擊事件后到2017年之前每年都會發(fā)生,那么A1處的撞擊則可能是后續(xù)5次異常事件的導(dǎo)火索.一種可能的推測是,這次撞擊位于尾部燃料箱附近,在每年的自轉(zhuǎn)加速段的幾個月時間里,在其中一個對稱性解中[4]太陽長時間照射尾部,持續(xù)的加熱引發(fā)受損區(qū)域的材料結(jié)構(gòu)變化從而造成燃料泄漏.而在2017年之后,由于殘留的燃料已完全釋放,再也沒有發(fā)生過轉(zhuǎn)速異常.我們可以預(yù)期在即將到來的任務(wù)之前不會再發(fā)生類似的加速事件.

4 最終的演化模型

除了A1外,在A4和A6也發(fā)生了特殊事件.圖2展示了北斗G2帆板參數(shù)的所有反演結(jié)果.在圖2(a)中,由S4和S5段反演得到的帆板正面總反射率相比之前陡然下降,而反面與正面的反射率比值(ratio of Back-to-Front,B/F ratio)與之前相比并沒有太大變化.這反映出目標所受到的力矩效應(yīng)減弱,但整體構(gòu)型沒有太大改變.一種可能的情況是在A4段目標的帆板受到了損壞.我們通過在S4和S5段中用單側(cè)帆板缺損25%的模型來反演,得到的正面總反射率與S1-S3的相當,且得到了與無損模型同樣的帆板角度和B/F ratio,驗證了我們的推測.這種損壞大概率來自目標快速旋轉(zhuǎn)造成的解體.在圖2(b)中也能看到該事件造成帆板角度分布的巨大間隔.此外,在A6段的起始處,即2016年6月29日(MJD=57568),NASA (National Aeronautics and Space Administration)公布稱北斗G2解體[9],并分裂出至少5塊(后續(xù)報道改為3 塊),但這些碎片至今沒有進TLE編目庫.圖2(a)中可以看出S6a和S6b的B/F ratio已經(jīng)發(fā)生了較大改變,說明衛(wèi)星確實發(fā)生了結(jié)構(gòu)性變化,比如其最大慣量主軸可能已不再沿著原體軸方向.理論上可以通過增加求解參數(shù)來確認損壞程度.但我們認為在圖1(e)-(f)中現(xiàn)有模型對S6a和S6b擬合后的外推結(jié)果與觀測數(shù)據(jù)都已有非常好的吻合度,新增參數(shù)可能會產(chǎn)生過擬合,所解出來結(jié)果的可靠性也不高.我們以S6b段反演得到的旋轉(zhuǎn)狀態(tài)為參考標準,將所建立的旋轉(zhuǎn)運動模型從S6a外推到S6b,并考慮了二者的不確定度,外推結(jié)果的標準差為轉(zhuǎn)軸赤經(jīng)0.5°,轉(zhuǎn)軸赤緯2.5°,轉(zhuǎn)速3.1×10-4s-1即0.11°·s-1.如果未來不再發(fā)生任何異常事件,該計算模型和計算參數(shù)能夠滿足對接時刻的旋轉(zhuǎn)狀態(tài)精度需求.最終,實踐21號順利完成了任務(wù).

圖2 北斗G2衛(wèi)星關(guān)于帆板的參數(shù)在不同段的反演結(jié)果.(a)帆板正面總反射率(Front total reflection coefficient,藍)和反面與正面反射率比值(ratio of back-to-front (B/F) reflection coefficients,紅).黑色為用S4和S5段在單側(cè)帆板缺損25%的模型來反演得到的正面總反射率.需要注意該反射率包含了模型求解中的多種因素,并不是真實反射系數(shù)[4].虛線標注了3次異常事件的時間.(b)兩個帆板角度δ1和δ2的解的分布.截取擬合殘差小于10-4 s-1 (S1-S5)和1.3×10-4 s-1 (S6a和S6b)的結(jié)果.Fig.2 Estimation results of parameters for the BeiDou-G2 satellite’s solar panels in all segments.(a) The front total reflection coefficient (blue) and the ratio of back-to-front (B/F) reflection coefficients (red).The black data represented the front total reflection coefficient estimated from the model with 25% of one solar panel being damaged in segments S4 and S5.The estimated reflection coefficient included various factors in the model simulation and was not the actual reflectivity or albedo[4].The dashed lines indicate the times of three abnormal events.(b) The distributions of panel angles δ1 and δ2 results.The fitting residuals of the results shown in this panel were less than 10-4 s-1 (for S1-S5) and 1.3×10-4 s-1 (for S6a and S6b).

致謝The author thanks Graz SLR station for providing photometry observation data in 2015.The author wish to thank Dr.Nilda Oklay,editor ofNature Communications,for her useful comments and efforts for this manuscript.感謝李龍華老師為本文賜名.感謝中科院DAN&NS小組對本文撰寫的幫助.謹以此文紀念愛子林澤瑜在本課題研究期間與我的相伴時光.