张 翔 HANCOCK Paul WAYTH Randall季江徽马月华

(1中国科学院紫金山天文台南京210008)

(2中国科学院行星科学重点实验室南京210008)

(3中国科学院大学北京100049)

(4 International Centre for Radio Astronomy Research,Curtin University,Perth WA6102)

1 引言

空间碎片又被称为太空垃圾,指的是分布在环绕地球轨道上(高度大约在100–40000 km之间)并已丧失功能的人造空间物体,包括废弃的卫星、火箭残骸等以及它们在碰撞、风化等过程中产生的碎片.空间碎片对航天器的飞行安全造成了重大威胁,碎片的撞击是航天器机械损伤的主要原因.例如2009年2月发生的俄罗斯废弃卫星Kosmos 2251与美国通信卫星Iridium 33相撞事故,导致Iridium 33报废,同时产生了大量新的空间碎片[1].因此,对空间碎片进行监测是十分必要的.

探测空间碎片的主要途径包括地基雷达和光学望远镜等.光学观测主要针对处在地球同步轨道(Geostationary Orbit,GEO,高度为35786 km)等高轨道的空间碎片,观测它们反射的太阳光,例如美国航空航天局(National Aeronautics and Space Administration,NASA)就用施密特望远镜对地球同步轨道上的空间碎片进行过巡天[2].雷达观测受距离的限制较大,但不受太阳光的影响,因而更适用于近地轨道(Low Earth Orbit,LEO,指2000 km高度以下的轨道范围)的空间碎片监测和成像,比如Goldstone雷达对部分3000 km以下空间碎片进行过观测[3],而Haystack雷达则监测过1000 km以下的空间环境[4−6].在例行的空间监测中,光学望远镜对地球同步轨道上空间碎片的探测极限大约为30 cm到1 m,而雷达对近地轨道空间碎片的探测极限大约为5–10 cm[7].

除了主动雷达之外,被动雷达(Passive Radar)用来观测空间碎片的可能性也在讨论中.被动雷达是双站雷达(Bistatic Radar)的一种特殊情况:雷达系统不包含发射器,通过接收目标散射的其他发射源的信号来进行探测.这些发射源可能是广播电台或者通讯台站等.比如Benson提出了利用全球卫星导航系统(代表为美国的全球定位系统,Global Positioning System,GPS)中的卫星作为发射站,并用地面天线接收空间碎片反射的卫星信号的方案[8].考虑到空间碎片散射的信号比较微弱,需要有大型射电天线参与地面信号接收[9].

近年来,Tingay等人提出了一种基于调频(Frequency Modulation,FM)广播信号的空间碎片探测方法[10]:使用默奇森大视场射电阵(Murchison Widefield Array,MWA,详见http://www.mwatelescope.org/)接收空间碎片散射的FM广播信号,进而实现对碎片的探测和追踪.在2014—2016年间,MWA进行了一系列针对流星自发辐射的观测[11],所使用波段与FM波段重叠.因此这部分数据也可以用来探测空间碎片.

本文第2节介绍MWA观测空间碎片的有利条件,包括MWA望远镜的一些参数及其探测空间碎片的能力估计;第3节介绍数据处理以及与空间天体相关的射电暂现源的筛选;第4节具体介绍观测到的射电暂现源;第5节对观测结果进行讨论;最后为小结.

2 观测设备

MWA是平方公里阵列[12](Square Kilometre Array,SKA)在低频波段的先导项目.它坐落于澳大利亚西部的默奇森射电天文台(Murchison Radio-astronomy Observatory,MRO),这里是平方公里阵列项目的两个台址之一(另一个是南非的北开普地区).默奇森射电天文台处于射电宁静区,受人类射频干扰(Radio Frequency Interference,RFI)的影响很小,尤其是在调频广播波段[13].

MWA由128个孔径阵(aperture array,见图1)组成,分布在3 km左右的范围内.它的等效面积大致为2500 m2,工作波段主要在80–300 MHz,带宽为30.72 MHz,最小时间分辨率为0.5 s,最小频率分辨率为10 kHz[14].对于空间碎片观测而言,在80 MHz左右,它有较大的视场范围(约2400 degree2)以及合适的角分辨率(约3′).

作为射电天线阵,MWA只能接收射电信号,没有发射功能.MWA附近也是射电宁静区,没有调频广播站.事实上,使用MWA作为被动雷达的探测方法正基于MWA距离地面广播站很遥远:如果MWA附近存在广播站,那么MWA直接从广播站接收的信号将比空间碎片反射的信号强得多,造成严重干扰.我们使用散布在澳大利亚各地(甚至其他国家和地区)的调频广播站作为发射器.澳大利亚西部的调频广播站见图2,其中调频广播站的地理位置数据来自澳大利亚通信和媒体管理局(Australian Communications and Media Authority,ACMA,详见https://www.acma.gov.au/).调频广播站工作波段大约为87.5–108 MHz.这些调频广播站的功率彼此间相差很大,从数瓦到数百千瓦不等.

图1 MWA核心区的数个孔径阵.图像来源:MWA/科廷大学Fig.1 Tiles from the core area of the MWA telescope.Image credit:MWA/Curtin

图2 澳大利亚西部的调频广播站分布.星号代表MWA所在地,圆点代表调频广播站所在地.纬度为南纬,经度为东经Fig.2 FM stations distributed in Western Australia.Location of the MWA is represented by a star,while the FM stations are represented by dots.Longitude is given as the east longitude,and latitude is give as the south latitude

2013年,Tingay等人曾经估算过MWA作为被动雷达对空间碎片的探测能力[10]:假设存在一块理想的空间碎片(完美导电的球体,半径为r),距离MWA为Rr,同时距离FM广播信号发射站Rt.FM广播站发射信号波长为λ.根据米氏散射公式[15],此空间碎片的雷达散射截面(Radar Cross Section,RCS,后面用σ来指代)大致可分为两种情况:(1)如果空间碎片尺寸远小于信号波长,即λ≫ 2πr,则雷达散射截面可以用瑞利散射来表示,σ=9πr2(2πr/λ)4;(2)如果空间碎片尺寸远大于信号波长,即λ≪ 2πr,则雷达散射截面可以用几何散射来表示,σ=πr2.考虑到小碎片的数量远多于大碎片,MWA接收到的某个碎片散射的信号大致可以由以下公式描述(根据Wills的双站雷达公式推导[16]):

其中,S是MWA接收到的谱流量密度(单位为Jy),Pt是FM广播站的有效发射功率(单位为kW),G为广播站发射在空间碎片方向的增益,B为发射信号带宽(单位为kHz),而ν为发射信号的中心频率(单位为MHz).基于(1)式,可以给出MWA对不同大小的碎片的大致观测距离范围(对1 m以上的碎片计算几何散射),见图3.根据进一步的模拟计算,MWA能够探测到1000 km高度处半径r>0.5 m的碎片.

图3 模拟空间碎片散射FM广播信号在MWA处产生的流量密度.假设空间碎片为理想导电球体,半径分别为0.3 m、1 m、3 m和10 m,距离MWA在400–2000 km 之间.图中红色虚线代表大致的MWA图像1σ噪声.Fig.3 Estimated flux density at the MWA as FM radio signals are scattered by space debris.Space debris are modeled as perfectly conducting spheres,with radii 0.3 m,1 m,3 m,and 10 m,at distances 400–2000 km from the MWA.The red dashed line gives approximately 1σ noise in the MWA images.

Palmer等人使用贝叶斯统计方法对MWA探测到国际空间站的数据进行了初轨计算,并将拟合出的轨道与国际空间站的两行轨道根数数据(Two Line Element,TLE)进行比较[17].在3 h的积分时间范围内,根据MWA拟合轨道推算的国际空间站位置与TLE根数推算的位置误差在数千米左右.

3 观测方法及数据处理

从2014年到2016年,MWA进行了一系列针对流星余迹辐射的观测[11].观测全部在夜间进行,总长度为322 h.考虑到我们不能预测流星何时出现在天空的什么位置,观测过程中MWA的所有天线指向天顶.观测频段为MWA的最低频段72–103 MHz,与澳大利亚的FM频段(87.5–108 MHz)部分重叠.根据观测模式的不同,MWA的时间分辨率从0.5 s到2 s不等,频率分辨率为40 kHz.

观测数据的处理分为3步:(1)使用Cotter[18]进行预处理.Cotter会将原始数据转化成通用天文软件包(Common Astronomy Software Applications,CASA)[19]的Measurement Set格式,同时Cotter可以去除数据中大部分的RFI.为了探测流星的自发辐射以及流星和空间碎片的反射,我们将数据分为两组,去除RFI的数据和保留RFI的数据.(2)使用明亮的射电点源对图像进行定标.常用的射电源包括Hydra A、Pictor A和3C 444.(3)使用WSClean[20]成像.为了提高目标信噪比,我们使用了8 s的积分时间,同时将前后两张图像相减得到差值图像.

我们使用一系列经验性步骤在未去除RFI的差值图像中寻找空间碎片反射调频广播信号产生的射电暂现事件.内容包括:(1)在图像中寻找移动暂现源形成的弧段.在8 s积分时间内,图像中心是固定在某个位置的,所以天空背景中大多数射电源在差值图像中会被消去.那些由于设备/空气扰动原因而未能完全消除的背景射电源也会显示成点源或者扩展源.但近地目标(包括流星、飞机、活动卫星、空间碎片等)相对MWA的方位变化较快,相当于移动的射电源,所以会在图像中产生弧段.(2)针对这些产生弧段的射电暂现事件,我们生成了2 s的积分图像,并寻找在连续4张图像中都出现的移动射电暂现源(见图4).这是为了排除流星:流星的电离余迹可以反射调频广播信号,但流星余迹形成所需的时间很短(大多数不到1 s).所以流星在MWA图像中可以形成弧段,但弧段会在形成的地方消散而不会移动.(3)对于前两步筛选剩下的暂现源,我们将它们和已知的近地轨道上的空间天体进行比较.我们使用北美防空司令部(North American Aerospace Defense Command,NORAD)提供的TLE轨道根数,并使用PyEphem软件包来计算相同时刻各个已知的卫星和空间碎片等(目前有TLE数据的空间天体数量大概在40000左右)相对MWA的位置,然后在MWA图像中标注出来(见图5),并与观测到的射电暂现源弧段进行比较.对于某个射电暂现源A和某个存在TLE数据的已知空间天体B,在连续4张2 s积分图像中,如果每张图的空间天体B的TLE轨道根数计算位置与暂现源A位置的角距离均在2◦以内,我们就认为暂现源A可能是由空间天体B反射FM广播信号而产生的.

图4 立方卫星DUCHIFAT 1的MWA差值观测图像.积分时间为2 s.图中黄色圆圈代表根据TLE数据计算出的在积分时间开始时刻DUCHIFAT 1的位置.图像中明亮弧段为MWA探测到的移动射电暂现源.Fig.4 MWA difference images of the CubeSat DUCHIFAT 1.Integration time is 2 s.Yellow circles in the image indicate the positions of DUCHIFAT 1 at the start of the integration time,as calculated with the TLE data.The bright curves in the image indicate a moving radio source observed by the MWA.

图4给出了一个MWA观测图像和轨道根数预报比较案例.图中较明亮的轨迹显示了某颗卫星或空间碎片相对天空背景的运动,轨迹中黑暗的部分是由前后两张图像相减造成的.图中黄色的圆圈则给出了在积分开始时间,根据TLE轨道根数估算的立方卫星DUCHIFAT 1所在位置.由此动态图像可以认为,MWA接收到了DUCHIFAT 1散射的FM广播信号.

4 观测结果

在322 h的观测数据中,我们一共确认了10起由近地空间天体引发的射电暂现事件.每起射电暂现事件持续数十秒到数分钟,信号源在图像中移动几度到几十度.与这些事件相关的近地空间天体大小从数百米到10 cm不等,高度则分布在500–2000 km之间.表1列出了这些射电暂现事件与相关天体,表中空间天体名称、近地点距离、远地点距离和雷达截面范围数据均来自NORAD.表中ISS代表国际空间站(International Space Station).雷达截面范围的定义如下:若RCS<0.1 m2,则雷达截面较小;若0.1 m2RCS1.0 m2,则雷达截面中等;若RCS>1.0 m2,则雷达截面较大.

根据引发射电暂现事件的近地空间天体的大小和距离,我们重点讨论以下3类事件.

(1)国际空间站及相关飞船,包括联盟号宇宙飞船(Soyuz),进步号宇宙飞船(Progress)和天鹅座宇宙飞船(Cygnus)等.国际空间站由于其尺寸远大于一般空间天体(超过100 m)且轨道高度较低(约400 km),是较为容易观测到反射FM信号的对象,在2013年测试MWA探测空间碎片能力的实验中就被观测到过[10].国际空间站的观测图像见图5.这类事件的特点是存在多条平行的运动轨迹线,可能与飞船有关,但也可能是空间站的结构造成的.

表1 空间天体引发的射电暂现事件.时间均为协调世界时.近地点和远地点距离的单位均为km.Table 1 Radio transient events caused by space objects.Date-Time in UTC;Apogee and perigee in km.

(2)ALOUETTE 2.它是一颗已经报废的加拿大卫星,大小约1 m.ALOUETTE 2射电暂现事件与其他事件的区别在于高度:被探测到反射FM信号时,它的高度约为2000 km.从图6可以看出,在积分时间相同(2 s)时,它的运动轨迹明显比空间站短.此外,由于轨道较高,ALOUETTE 2事件持续时间(4 min)比其他事件(约数十秒)要长,空间天体可观测范围(数十度)也比其他天体(数度)要大.

从ALOUETTE 2事件中可以看出,使用差值观测图像寻找空间碎片的方法可能存在高度限制.我们使用差值图像来降低背景噪声(相邻积分时间图像相减),但在空间天体比较高的情况下,若积分时间较短,则天体的移动不明显.这可以导致空间天体在差值图像中被消去.换言之,使用反射FM信号探测空间碎片,需要根据目标碎片的高度不同而设置不同的积分时间.

(3)DUCHIFAT 1和UKUBE 1.它们都是立方卫星,其中DUCHIFAT 1的尺寸只有0.1 m.这表明MWA拥有探测到米级以下空间碎片的潜力.DUCHIFAT 1的观测图像已在图4中给出.

图6 ALOUETTE 2的MWA差值观测图像.积分时间为2 s.图中黄色圆圈代表根据TLE数据计算出的在2 s开始时刻各空间天体的位置.Fig.6 MWA difference image of ALOUETTE 2.Integration time is 2 s.Yellow circles in the image indicate the positions of the space objects at the start of the integration time,as calculated with the TLE data.

5 讨论

与其他观测方法相比,使用调频广播信号探测空间碎片存在分辨率较低的问题.由于频率较低而波长较长(在100 MHz处波长为3 m),MWA作为被动雷达观测空间碎片的空间分辨率大约为4.5′(在600 km高度处相当于0.78 km),导致空间碎片位置数据中出现沿迹误差(along-track error)和垂迹误差(cross-track error).除此之外,为了提高观测图像信噪比,我们可能需要对观测数据在时间上进行积分,导致时间分辨率降低.考虑到空间碎片在观测图像上的移动,这主要会引入沿迹误差.

我们用两个事例—DUCHIFAT 1和ALOUETTE 2来讨论MWA作为被动雷达探测空间碎片的误差大小.在计算中,我们将MWA观测到的空间天体位置与使用TLE轨道根数计算出的空间天体位置相比较,并假设使用TLE轨道根数计算出的空间天体位置与真实值的误差远小于MWA观测得到的位置与真实值的误差.

图7中给出了DUCHIFAT 1的观测结果和误差估计.上方图中圆点代表在10 s的连续观测中,每2 s积分时间开始时DUCHIFAT 1在MWA观测图像中的位置,而误差范围是由MWA的空间分辨率和积分时间来估算的.图中红色三角形给出的是在每2 s积分开始时刻根据TLE根数推算的DUCHIFAT 1在MWA观测图像中的位置.由图7可以看出,对于600 km高度处的DUCHIFAT 1而言,MWA空间和时间分辨率导致的方位误差大概在1.5◦左右,而TLE根数推算出的位置与此相符.在已知空间天体在MWA图像上的位置和TLE根数计算的位置之间的角距离和空间天体高度的前提下,可以估算MWA作为被动雷达观测到的空间天体位置与其TLE根数推算得到的位置之间的距离,见图7中下方图.但此处计算的距离没有考虑视向方向上的误差,所以会比真实值小.从图中可以看出,对DUCHIFAT 1而言,忽略视向方向的MWA观测位置和TLE根数推算位置的差距在10 km以内.

图7 从MWA的方位观测DUCHIFAT 1的运动轨迹.上图:圆点代表MWA图像中与DUCHIFAT 1对应的射电暂现源在10 s观测时间内的位置,误差范围根据MWA的空间和时间分辨率给出.采样时间为2 s.三角形给出的是在相同时间根据TLE轨道根数推算的DUCHIFAT 1在MWA图像中的位置.下图:在10 s观测时间中,估算的MWA观测到的空间天体位置与根据TLE根数推算得到的位置之间的距离.Fig.7 The path of DUCHIFAT 1 as seen by the MWA.Top:the dots indicate the positions of radio transients caused by DUCHIFAT 1 during 10 s in the MWA images.Errors are given by estimation with temporal and spatial resolution of the MWA.Sampling timescale is 2 s.The triangles represent the positions of DUCHIFAT 1 in the MWA images calculated with the TLE data.Bottom:estimated distance between the positions of DUCHIFAT 1 as observed by the MWA and calculated with the TLE data,during 10 s observation.

ALOUETTE 2事件则有所不同.在图8中,我们给出了ALOUETTE 2在14 s时间中的观测结果和误差估计.由于在观测时刻ALOUETTE 2所在高度较高(约2000 km),运动的角速度较小,导致积分时间引入的沿迹误差比DUCHIFAT 1小一些.我们估算MWA观测2000 km高度的空间碎片的方位误差大致小于1◦.从图8可以看出:根据TLE轨道根数估算的ALOUETTE 2位置在MWA观测的大致误差范围之外,相差0.1◦左右.这可能与TLE数据自身的误差,使用的Pyephem软件包在根据TLE数据计算空间天体的方位过程中产生的误差以及我们使用的ALOUETTE 2的TLE轨道根数历元与观测时间相差较大有关.ALOUETTE 2的轨道根数历元与观测时间相差11 h.在DUCHIFAT 1事件中这一差距为4 h.

图8 从MWA的方位观测ALOUETTE 2的运动轨迹.上图:圆点代表MWA图像中与ALOUETTE 2对应的射电暂现源在14 s观测时间内的位置,误差范围根据MWA的空间和时间分辨率给出.采样时间为2 s.三角形给出的是在相同时间根据TLE轨道根数推算的ALOUETTE 2在MWA图像中的位置.下图:在14 s观测时间中,估算的MWA观测到的空间天体位置与根据TLE根数推算得到的位置之间的距离.Fig.8 The path of ALOUETTE 2 as seen by the MWA.Top:the dots indicate the positions of radio transients caused by ALOUETTE 2 during 14 s in the MWA images.Errors are given by estimation with temporal and spatial resolution of the MWA.Sampling timescale is 2 s.The triangles represent the positions of ALOUETTE 2 in the MWA images calculated with the TLE data.Bottom:estimated distance between the positions of ALOUETTE 2 as observed by the MWA and calculated with the TLE data,during 14 s observation.

此外还值得注意的一点:在我们手动寻找到的10次射电暂现现象中,DUCHIFAT 1占了4次,国际空间站和飞船占了3次.空间站被重复观测到可能和它的尺寸有较大关系;而DUCHIFAT 1作为立方卫星,被多次观测到可能有其他原因.为此我们统计了FM信号反射射电暂现事件的空间分布(天顶距和方位角),结果见图9.从图9中可以看出,有3次DUCHIFAT 1事件的方位角非常接近.我们猜测一种可能的解释是反射的FM信号来自同一个FM广播站,也就是说DUCHIFAT 1被多次观测到是因为它的位置和FM广播站以及MWA形成了合适的几何结构.这个猜测是根据流星雨观测的结果类比出的:在流星雨的被动雷达观测中,如果已知发射机和接收机的位置,则可以估算流星能够被观测到的天空范围[21].换言之,发射机、接收机以及流星三者的相互位置关系会影响流星被探测到的几率.在这些DUCHIFAT 1事件中,如果在MWA某个方向的合适位置上存在一个FM广播站,或许可以解释DUCHIFAT 1在这个方向被多次探测到.

图9 MWA观测到的反射FM广播信号射电暂现事件的空间分布.空间位置由天顶距和方位角给出.Fig.9 Space distribution of radio transient events caused by the reflected FM signal,as observed by the MWA.Positions of the events are given in Azimuth and Zenith Angle.

6 总结与展望

我们使用了默奇森大视场射电阵搜寻了空间天体反射的FM广播信号,探测到10起反射FM信号引起的射电暂现事件.将这些事件与轨道根数计算得出的空间天体位置相对照,我们确认这些事件是由数颗卫星引起的.这些事件证明了MWA拥有探测600 km高度处0.1 m量级尺寸空间碎片的潜力,而观测的位置误差大概在10 km左右.

本文中的数据处理基于MWA流星观测项目,因此许多参数的选择相对空间碎片观测并不是最优的.例如积分时间可以根据目标碎片高度而修改,而观测带宽可以由MWA周围较强的FM广播波段确定.此外,本文中介绍的所有射电暂现事件都是手动选出的,尝试使用图像搜索软件也许能得到更多的事件.

根据估算,本文中MWA观测空间天体的位置误差大部分是由数据处理中采用了较长的积分时间导致的,在进一步的数据处理中可以考虑减小积分时间,甚至不对时间做积分.然而,积分时间减小会引起信噪比的降低,这方面的平衡也是我们需要考虑的.

在本文中,MWA作为被动雷达只测量了空间天体的方位,没有测距和测速数据.这是因为对某起空间天体导致的射电暂现事件而言,我们尚未确认发射调频广播信号的具体广播站.我们下一步可以根据观测到的射电暂现源的频率分布和当地广播站位置分布得到具体的广播站,进而将MWA接收到的信号与发射信号相比较,提高观测的精度.

探测近地轨道空间碎片的主要方式是主动雷达,探测极限在5–10 cm左右.到目前为止,MWA作为被动雷达展现出来的探测能力尚未超过这一极限.但使用低频射电望远镜作为被动雷达的空间碎片探测方法,与建设主动雷达站相比,是一种成本较低的探测方式.另外,随着SKA项目的推进,默奇森天文台正在建造新的SKA-Low天线阵.此天线阵同样覆盖FM观测波段,同时灵敏度相对MWA会大幅提高.在可预期的未来,SKA-Low是进行空间碎片监测的合适工具.