关于射电望远镜台址航空信号处理方法的研究∗

许红瑞陈卯蒸 刘 奇 尹 航 王 军刘 璇袁 力

(1中国科学院新疆天文台乌鲁木齐830011)

(2中国科学院射电天文重点实验室南京210008)

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

1 引言

H I(0.235–0.386 z)和1665 OH(0.448–0.624 z)是研究天体移动及规律预测的重要谱线,但由于红移,谱线频率与航空信号频率发生重叠,如表1所示.民用航空信号(civil aviation signal,以下简称CAS)的工作波段为1–1.5 GHz之间约130 MHz带宽的L波段.在正常的射电天文观测中,CAS信号会从天线旁瓣进入接收机系统,降低系统信噪比和观测数据的有效性,增加数据处理的复杂性,严重影响了射电天文在此频段内的观测效果,降低大型射电科学仪器的使用价值[1−2].鉴于民航航线的不断增多和原有航线的调整,在新一代射电望远镜的选址工作中,虽然考虑了对地面无线电信号的电磁兼容和电磁干扰,但是对来自CAS的干扰却难以避免.

CAS主要来自民用航空飞机的机载广播式自动相关监视系统(Automatic Dependent Surveillance-Broadcast,以下简称ADS-B)和机载测距机系统(Distance Measuring Equipment,以下简称DME),ADS-B和DME信号占用L波段带宽如表1所示.ADS-B的作用距离可达300 km以上,系统发射功率为51–57 dBm,信号编码方式为脉冲编码,比特率为1 Mb/s[2].DME系统机载询问器发射频率在1025–1150 MHz,带宽占用125 MHz,信号带宽1–3 MHz,响应时间128µs[3].即ADS-B和机载DME产生的电磁干扰信号具有瞬态、极化和高辐射强度等特征.目前射电观测数据的消干扰方法为屏蔽子积分频谱中有干扰的通道,基于现有的消干扰方法,在干扰信号越多的情况下,有用的频谱通道越少,获得的天文数据的有用信息越少[4].国内外针对由瞬时脉冲信号造成的射频干扰多集中在后期处理,Willem Baan团队提出避免干扰信号比抑制干扰信号更有意义,通过利用自然地形,提前获知所观测区域有哪些干扰频谱存在,尽量避免或者消除干扰,从而有效建立观测宁静区[5].

为避免射电望远镜观测H I(0.235–0.386 z)和1665 OH(0.448–0.624 z)[6−8]红移过程中出现的CAS,本文利用已开发的软件无线电平台,通过对ADS-B报文解析获得航迹样本点.基于航迹样本点分布的时间特点,采用最小二乘多项式拟合,获得航迹样本点在时间上的分布高峰时段,结合获得的高峰时段对航迹样本点进行基于球面距离的二分k-均值聚类算法分析,从而获得射电天文台址周边高峰时段飞机分布的主要区域.为安排射电天文观测、避开CAS、建设射电天文观测宁静区及消干扰策略提供依据.

表1 飞机导航信号及红移谱线在1025–1150 MHz频段占用带宽Table 1 Frequency occupation of the aircraft navigation signal and the redshift spectrum in 1025–1150 MHz

2 数据存储格式及解析过程

本文基于已开发搭建的软件无线电平台,以新疆天文台本部大楼为测试点,经过测试解析,共获得航迹信息12万条,数据存储格式如表2所示,其中第1列是飞机唯一的ICAO(国际民用航空组织)地址,Seen表示的是接收到飞机信号的实时时间,Altitude、Latitude和Longitude分别代表飞机此时的海拔高度、纬度和经度.ADS-B信号经天线接收,经过射频处理之后成为基带信号,经传输到达PC端进行报文解析,获得航迹信息包括接收时间、地址码、航班号、经纬度、高度等多个报文项.解析流程如图1所示,其中,RF为射频前端,AD为数模转换子板,FPGA为现场可编程阵列,USRP是通用软件无线电外设.

表2 ADS-B报文解析存储格式Table 2 ADS-B message parse storage format

图1 ADS-B报文接收解析框图Fig.1 Receiving and parse block diagram of ADS-B message

3 飞机到测站距离及到达测站的干扰强度估算

3.1 计算空间距离

由于飞机在飞行过程中是动态的,接收到的点均为离散点,为了更直观地表现飞机到测试点距离在时间上的分布,计算C、D两点距离,其中C为飞机所处位置,其纬度、经度、海拔分别为a1、c1、b1,D为测试点所处位置,其纬度、经度、海拔分别为a2、c2、b2.利用(1)–(3)式,分别将纬度、经度、海拔转换为3维空间坐标x1、y1、z1和x2、y2、z2,利用(4)式计算C、D两点距离d,根据距离画出飞机在时间上的分布.

3.2 估算干扰强度

本文采用自由空间损耗模型估算信号强度,空间损耗公式为:

其中f为频率,单位为MHz,d为距离,单位为km.

3.3 飞机到测站距离随着时间分布及到测站信号强度估算结果

根据(4)式计算了一天中飞机到测站距离随着时间的分布图,如图2所示.由图2可知:飞机到测站距离分布集中在20 km左右,分布范围在10–80 km.相应地,若飞机在视线范围,到达测站的自由空间损耗在113–131 dB之间,若因地面等遮挡,飞机不在视线范围内,损耗一般会比同距离自由空间更大些,由此可以估算航空信号在观测站的干扰强度.

图2 飞机到测站距离随时间分布图Fig.2 The distribution of the distance between plane and station with time

4 基于最小二乘法和二分k-均值聚类算法的航迹分析方法

航线的主要分布时段和区域是众多航线在时间空间上反映出来的主要特征,提取航线高峰时段分布区域的技术路线如图3所示.获取测试点航线高峰时段主要区域基本思路为:首先根据获得的航迹点数量和时间的关系,分析每天接收的航迹点是否稳定,针对航迹点统计量在24 h的分布情况,采用最小二乘法的多项式拟合获得每天的高峰时段.在获得高峰时段后,航迹点按照高峰时段筛选出经纬度信息,之后应用基于球面距离的二分k-均值聚类算法对高峰时段的经纬度信息进行聚类,获得主要航线分布范围.

4.1 高峰时段分析

4.1.1 拟合原理

在2017年3月5日到2017年3月17日之间,每天接收解析信息总次数如表3所示.在不考虑平台性能和天气原因的情况下,总的来看,每天接收信号数量大致是稳定的.在此基础上,分析一天的高峰时段.本文认为每天的高峰时段应该是相同的,通过统计3月5日到3月17日之间每天24 h中每小时接收信号的总频次,得到13组数据点,每组24个统计点.首先计算对应各统计点每天平均样本数,得到一组24个平均数据,再按这组平均数据进行基于最小二乘法的多项式拟合,得到航迹样本点随时间的分布趋势.

图3 技术路线图Fig.3 Technology roadmap

表3 13组样本的统计量Table 3 Statistics of 13 samples

4.1.2 最小二乘多项式拟合

给定离散数据点(xi,yi)(i=0,1,2,···,m),xi代表时间点,yi表示对应时间点接收的统计量,Φ为次数不超过n(n≤m)的多项式构成的函数类,通过给定的离散数据点建立拟合函数为:使得(6)式的平方和I最小.

满足(6)式的pn(x)为最小二乘拟合多项式.

即

(8)式关于a0,a1,···,an的线性方程组,用矩阵表示为:

方程组(9)的矩阵是一个对称正定矩阵,故存在唯一解.从(9)式中解出ak(k=0,1,···,n),从而可得多项式:

4.1.3 相关系数及其显著性检验

将平均样本点(xi,yi)做多项式拟合时,还并不太了解x与y之间关系的密切程度.为此要用相关系数r进行判断,其定义为:

4.2 基于球面距离的二分k-均值航迹聚类算法

本文目的是将目标航迹进行聚类[9]分析,将相似航迹点聚在一起,在众多杂乱无章的航迹中找到出现航迹概率较大的区域范围.虽然k均值方法也可以很好地将数据进行聚类分析,但是二分k-均值聚类算法是对k均值聚类算法的改进,一方面可以克服k均值聚类算法收敛于局部最小值的问题,另一方面可以加快运算速度.因本文数据为经纬度信息,为了保证方法的普遍适用性,本文采用基于球面距离的二分k-均值聚类算法对航迹样本点进行处理.

4.2.1 计算球面距离

由于我们的数据主要是经纬度,但这些信息对于距离计算并不足够.在北极附近每走几米的经度变化可能达到10◦;而在赤道附近走相同的距离,带来的变化可能仅零点几米.为了统一规范化距离的计算,此处采用球面余弦定理来计算两点之间的距离.A点纬度、经度为(β1,α1),B点纬度、经度为(β2,α2),R为地球半径,则A、B两点间的球面距离d1为:

4.2.2 航迹聚类处理算法

本文提出了基于聚类的航迹分析算法,其输入包括,输入1:航迹集合Q={q1,q2,q3,···,qn},其中qi(i=1,2,···,n)代表航迹点;输入2: 指定k为希望得到的簇数目.输出为航迹聚类结果,算法具体执行步骤描述如下:利用二分k-均值算法对航迹进行分析,在算法中使用的距离为(12)式中的球面距离,得到航迹聚类结果,其中二分k-均值算法伪代码描述为:

(1)将所有点看成一个簇;

(2)当簇数目小于k时,对于每一个簇,进行:

(i)计算总误差;

(ii)在给定的簇上面进行k均值聚类;

(iii)计算将该簇一分为二之后的总误差;

(3)选择误差最小的那个簇进行划分操作.

5 实验结果分析与方法验证

5.1 基于最小二乘法的曲线拟合结果及分析

基于2017年3月5日到2017年3月17日之间一天24 h的13组统计数据(每组数据有24个统计点),按一天24个统计点,先计算对应各统计点每天平均样本数,再按这组平均数据进行基于最小二乘法的多项式拟合获得的趋势走向如图4所示.由图中连续的拟合曲线可以看出,一天之中的高峰时段分布在UTC时间的3点、9点和15点内,图中离散点为2017年3月5日到2017年3月17日之间,24个统计点所对应的实测统计值.

经计算对应时间点的平均样本数的相关性,以r=0.5为阈值,其中相关系数在阈值以上的占92%,因此我们认为得到的高峰时段的数据是可靠的,基于此,进行航迹聚类的实验.

图4 2017年3月5日—2017年3月17日之间每日平均样本数拟合分布趋势Fig.4 Fitting of the daily average sample counting distribution between 2017-03-05—2017-03-17

5.2 高峰时段航迹分布范围聚类划分

采用实测数据对航迹聚类算法进行测试,数据分为6组,前3组作为测试划分范围,后3组验证方法的有效性.第1组为3月5日–3月11日UTC时间3点内所接收的全部航迹点,第2组为3月5日–3月11日UTC时间9点内所接收的全部航迹点,第3组为3月5日–3月11日UTC时间15点内所接收的全部航迹点,第4组为3月12日–3月17日UTC时间3点内所接收的全部航迹点,第5组为3月12日–3月17日UTC时间9点内所接收的全部航迹点,第6组为3月12日–3月17日UTC时间15点内所接收的全部航迹点,由1、2、3组航迹数据聚类结果分别如图5(左侧)、6(左侧)、7(左侧)所示.根据1、2、3组数据航迹聚类结果我们可以划定高峰时段航迹出现的主要范围区域.图中红色圆点代表新疆天文台总部大楼,蓝色圆点代表航迹点,紫色、红色、绿色圈定的区域为在此时段内航迹点主要聚集范围.

图5 UTC 3点航迹聚类(左)和验证航迹聚类(右)Fig.5 Track clustering during UTC 3(left)and track clustering verifying during UTC 3(right)

5.3 方法验证

根据实验划定高峰时段范围区域,将后3组高峰时段航迹点作为验证数据,聚类结果如图5(右侧)、6(右侧)、7(右侧)所示.由图中可以看出,后3组航迹点经聚类后多数落入划定区域范围内.统计后3组高峰时段航迹点出现在划定范围区域的概率,我们将高峰时段接收的总的航迹点的集合称作S={s1,s2···sn},落在区域范围外的航迹点的集合称为落在区域范围内概率公式为:验证测试所得结果如表4.经过聚类算法的处理结果表明:航迹在高峰时段的范围比较确定,符合实际航迹分布特点.

图6 UTC 9点航迹聚类(左)和验证航迹聚类(右)Fig.6 Track clustering during UTC 9(left)and track clustering verifying during UTC 9(right)

图7 UTC 15点航迹聚类(左)和验证航迹聚类(右)Fig.7 Track clustering during UTC 15(left)and track clustering verifying during UTC 15(right)

表4 概率统计Table 4 Probability statistics

6 结论

以新疆天文台本部大楼为测试点,针对在测试点接收解析的时间、经纬度的航迹样本数据,对数据按照时间节点初步统计筛选后,根据数据分布的时间特点,选用基于最小二乘法的多项式对数据进行拟合,获得航迹样本在一天的分布趋势,经计算相关系数,验证了用最小二乘法的多项式分析趋势的可靠性.

根据高峰时段的分布,将获得的航迹样本点分为两组,一组用于测试,另一组用于验证方法的有效性.在对第一组高峰时段数据进行处理之后获得高峰时段航迹样本主要分布范围后,利用相同的方法对第二组数据进行航迹样本的聚类处理,经验证航迹样本点高达96%落在划定的区域范围内,验证了该方法的有效性及实用性.通过数据拟合可有效地分析一天之中的高峰时段的分布,根据拟合结果对高峰时段航迹进行处理,划定航迹出现的主要范围区域,即可获得一天之后高峰时段的航迹分布范围,从而为射电天文观测时段提供有效的参考,最大可能地避开射电天文观测中CAS对射电天文观测的影响和对射电望远镜的危害,为规划建设中的射电天文台址划分宁静区提供与相关部门交涉的依据.

[1]苗可可,王壮,程翥,等.天文研究与技术,2015,12:433

[2]Jones S R.Air Traffic Control Quaterly,2003,11:225

[3]Robert J F,Zhang Q,Zheng Y,et al.ApJ,2005,129:2940

[4]张群涛,刘奇,孙正文,等.科学技术与工程,2016,1:191

[5]李耀华.基于LMS的射电天文大天线抗干扰技术的研究与实现.长沙:国防科技大学,2011:8

[6]吴忠祖,Haynes M P,Giovanlli R,等.天文学报,2015,56:112

[7]Chang T C,Pen U L,Bandura K,et al.Nature,2010,466:463

[8]Steven W E,Grant H.ApJS,2003,147:167

[9]Maxime G,Ashok N S.IEEE Transactions on Intelligent Transportation Systems,2011,12:1511