許紅瑞陳卯蒸 劉 奇 尹 航 王 軍劉 璇袁 力

(1中國科學(xué)院新疆天文臺烏魯木齊830011)

(2中國科學(xué)院射電天文重點實驗室南京210008)

(3中國科學(xué)院大學(xué)北京100049)

1 引言

H I(0.235–0.386 z)和1665 OH(0.448–0.624 z)是研究天體移動及規(guī)律預(yù)測的重要譜線,但由于紅移,譜線頻率與航空信號頻率發(fā)生重疊,如表1所示.民用航空信號(civil aviation signal,以下簡稱CAS)的工作波段為1–1.5 GHz之間約130 MHz帶寬的L波段.在正常的射電天文觀測中,CAS信號會從天線旁瓣進(jìn)入接收機(jī)系統(tǒng),降低系統(tǒng)信噪比和觀測數(shù)據(jù)的有效性,增加數(shù)據(jù)處理的復(fù)雜性,嚴(yán)重影響了射電天文在此頻段內(nèi)的觀測效果,降低大型射電科學(xué)儀器的使用價值[1?2].鑒于民航航線的不斷增多和原有航線的調(diào)整,在新一代射電望遠(yuǎn)鏡的選址工作中,雖然考慮了對地面無線電信號的電磁兼容和電磁干擾,但是對來自CAS的干擾卻難以避免.

CAS主要來自民用航空飛機(jī)的機(jī)載廣播式自動相關(guān)監(jiān)視系統(tǒng)(Automatic Dependent Surveillance-Broadcast,以下簡稱ADS-B)和機(jī)載測距機(jī)系統(tǒng)(Distance Measuring Equipment,以下簡稱DME),ADS-B和DME信號占用L波段帶寬如表1所示.ADS-B的作用距離可達(dá)300 km以上,系統(tǒng)發(fā)射功率為51–57 dBm,信號編碼方式為脈沖編碼,比特率為1 Mb/s[2].DME系統(tǒng)機(jī)載詢問器發(fā)射頻率在1025–1150 MHz,帶寬占用125 MHz,信號帶寬1–3 MHz,響應(yīng)時間128μs[3].即ADS-B和機(jī)載DME產(chǎn)生的電磁干擾信號具有瞬態(tài)、極化和高輻射強(qiáng)度等特征.目前射電觀測數(shù)據(jù)的消干擾方法為屏蔽子積分頻譜中有干擾的通道,基于現(xiàn)有的消干擾方法,在干擾信號越多的情況下,有用的頻譜通道越少,獲得的天文數(shù)據(jù)的有用信息越少[4].國內(nèi)外針對由瞬時脈沖信號造成的射頻干擾多集中在后期處理,Willem Baan團(tuán)隊提出避免干擾信號比抑制干擾信號更有意義,通過利用自然地形,提前獲知所觀測區(qū)域有哪些干擾頻譜存在,盡量避免或者消除干擾,從而有效建立觀測寧靜區(qū)[5].

為避免射電望遠(yuǎn)鏡觀測H I(0.235–0.386 z)和1665 OH(0.448–0.624 z)[6?8]紅移過程中出現(xiàn)的CAS,本文利用已開發(fā)的軟件無線電平臺,通過對ADS-B報文解析獲得航跡樣本點.基于航跡樣本點分布的時間特點,采用最小二乘多項式擬合,獲得航跡樣本點在時間上的分布高峰時段,結(jié)合獲得的高峰時段對航跡樣本點進(jìn)行基于球面距離的二分k-均值聚類算法分析,從而獲得射電天文臺址周邊高峰時段飛機(jī)分布的主要區(qū)域.為安排射電天文觀測、避開CAS、建設(shè)射電天文觀測寧靜區(qū)及消干擾策略提供依據(jù).

表1 飛機(jī)導(dǎo)航信號及紅移譜線在1025–1150 MHz頻段占用帶寬Table 1 Frequency occupation of the aircraft navigation signal and the redshift spectrum in 1025–1150 MHz

2 數(shù)據(jù)存儲格式及解析過程

本文基于已開發(fā)搭建的軟件無線電平臺,以新疆天文臺本部大樓為測試點,經(jīng)過測試解析,共獲得航跡信息12萬條,數(shù)據(jù)存儲格式如表2所示,其中第1列是飛機(jī)唯一的ICAO(國際民用航空組織)地址,Seen表示的是接收到飛機(jī)信號的實時時間,Altitude、Latitude和Longitude分別代表飛機(jī)此時的海拔高度、緯度和經(jīng)度.ADS-B信號經(jīng)天線接收,經(jīng)過射頻處理之后成為基帶信號,經(jīng)傳輸?shù)竭_(dá)PC端進(jìn)行報文解析,獲得航跡信息包括接收時間、地址碼、航班號、經(jīng)緯度、高度等多個報文項.解析流程如圖1所示,其中,RF為射頻前端,AD為數(shù)模轉(zhuǎn)換子板,FPGA為現(xiàn)場可編程陣列,USRP是通用軟件無線電外設(shè).

表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 飛機(jī)到測站距離及到達(dá)測站的干擾強(qiáng)度估算

3.1 計算空間距離

由于飛機(jī)在飛行過程中是動態(tài)的,接收到的點均為離散點,為了更直觀地表現(xiàn)飛機(jī)到測試點距離在時間上的分布,計算C、D兩點距離,其中C為飛機(jī)所處位置,其緯度、經(jīng)度、海拔分別為a1、c1、b1,D為測試點所處位置,其緯度、經(jīng)度、海拔分別為a2、c2、b2.利用(1)–(3)式,分別將緯度、經(jīng)度、海拔轉(zhuǎn)換為3維空間坐標(biāo)x1、y1、z1和x2、y2、z2,利用(4)式計算C、D兩點距離d,根據(jù)距離畫出飛機(jī)在時間上的分布.

3.2 估算干擾強(qiáng)度

本文采用自由空間損耗模型估算信號強(qiáng)度,空間損耗公式為:

其中f為頻率,單位為MHz,d為距離,單位為km.

3.3 飛機(jī)到測站距離隨著時間分布及到測站信號強(qiáng)度估算結(jié)果

根據(jù)(4)式計算了一天中飛機(jī)到測站距離隨著時間的分布圖,如圖2所示.由圖2可知:飛機(jī)到測站距離分布集中在20 km左右,分布范圍在10–80 km.相應(yīng)地,若飛機(jī)在視線范圍,到達(dá)測站的自由空間損耗在113–131 dB之間,若因地面等遮擋,飛機(jī)不在視線范圍內(nèi),損耗一般會比同距離自由空間更大些,由此可以估算航空信號在觀測站的干擾強(qiáng)度.

圖2 飛機(jī)到測站距離隨時間分布圖Fig.2 The distribution of the distance between plane and station with time

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

航線的主要分布時段和區(qū)域是眾多航線在時間空間上反映出來的主要特征,提取航線高峰時段分布區(qū)域的技術(shù)路線如圖3所示.獲取測試點航線高峰時段主要區(qū)域基本思路為:首先根據(jù)獲得的航跡點數(shù)量和時間的關(guān)系,分析每天接收的航跡點是否穩(wěn)定,針對航跡點統(tǒng)計量在24 h的分布情況,采用最小二乘法的多項式擬合獲得每天的高峰時段.在獲得高峰時段后,航跡點按照高峰時段篩選出經(jīng)緯度信息,之后應(yīng)用基于球面距離的二分k-均值聚類算法對高峰時段的經(jīng)緯度信息進(jìn)行聚類,獲得主要航線分布范圍.

4.1 高峰時段分析

4.1.1 擬合原理

在2017年3月5日到2017年3月17日之間,每天接收解析信息總次數(shù)如表3所示.在不考慮平臺性能和天氣原因的情況下,總的來看,每天接收信號數(shù)量大致是穩(wěn)定的.在此基礎(chǔ)上,分析一天的高峰時段.本文認(rèn)為每天的高峰時段應(yīng)該是相同的,通過統(tǒng)計3月5日到3月17日之間每天24 h中每小時接收信號的總頻次,得到13組數(shù)據(jù)點,每組24個統(tǒng)計點.首先計算對應(yīng)各統(tǒng)計點每天平均樣本數(shù),得到一組24個平均數(shù)據(jù),再按這組平均數(shù)據(jù)進(jìn)行基于最小二乘法的多項式擬合,得到航跡樣本點隨時間的分布趨勢.

圖3 技術(shù)路線圖Fig.3 Technology roadmap

表3 13組樣本的統(tǒng)計量Table 3 Statistics of 13 samples

4.1.2 最小二乘多項式擬合

給定離散數(shù)據(jù)點(xi,yi)(i=0,1,2,···,m),xi代表時間點,yi表示對應(yīng)時間點接收的統(tǒng)計量,Φ為次數(shù)不超過n(n≤m)的多項式構(gòu)成的函數(shù)類,通過給定的離散數(shù)據(jù)點建立擬合函數(shù)為:使得(6)式的平方和I最小.

滿足(6)式的pn(x)為最小二乘擬合多項式.

即

(8)式關(guān)于a0,a1,···,an的線性方程組,用矩陣表示為:

方程組(9)的矩陣是一個對稱正定矩陣,故存在唯一解.從(9)式中解出ak(k=0,1,···,n),從而可得多項式:

4.1.3 相關(guān)系數(shù)及其顯著性檢驗

將平均樣本點(xi,yi)做多項式擬合時,還并不太了解x與y之間關(guān)系的密切程度.為此要用相關(guān)系數(shù)r進(jìn)行判斷,其定義為:

4.2 基于球面距離的二分k-均值航跡聚類算法

本文目的是將目標(biāo)航跡進(jìn)行聚類[9]分析,將相似航跡點聚在一起,在眾多雜亂無章的航跡中找到出現(xiàn)航跡概率較大的區(qū)域范圍.雖然k均值方法也可以很好地將數(shù)據(jù)進(jìn)行聚類分析,但是二分k-均值聚類算法是對k均值聚類算法的改進(jìn),一方面可以克服k均值聚類算法收斂于局部最小值的問題,另一方面可以加快運(yùn)算速度.因本文數(shù)據(jù)為經(jīng)緯度信息,為了保證方法的普遍適用性,本文采用基于球面距離的二分k-均值聚類算法對航跡樣本點進(jìn)行處理.

4.2.1 計算球面距離

由于我們的數(shù)據(jù)主要是經(jīng)緯度,但這些信息對于距離計算并不足夠.在北極附近每走幾米的經(jīng)度變化可能達(dá)到10?;而在赤道附近走相同的距離,帶來的變化可能僅零點幾米.為了統(tǒng)一規(guī)范化距離的計算,此處采用球面余弦定理來計算兩點之間的距離.A點緯度、經(jīng)度為(β1,α1),B點緯度、經(jīng)度為(β2,α2),R為地球半徑,則A、B兩點間的球面距離d1為:

4.2.2 航跡聚類處理算法

本文提出了基于聚類的航跡分析算法,其輸入包括,輸入1:航跡集合Q={q1,q2,q3,···,qn},其中qi(i=1,2,···,n)代表航跡點;輸入2: 指定k為希望得到的簇數(shù)目.輸出為航跡聚類結(jié)果,算法具體執(zhí)行步驟描述如下:利用二分k-均值算法對航跡進(jìn)行分析,在算法中使用的距離為(12)式中的球面距離,得到航跡聚類結(jié)果,其中二分k-均值算法偽代碼描述為:

(1)將所有點看成一個簇;

(2)當(dāng)簇數(shù)目小于k時,對于每一個簇,進(jìn)行:

(i)計算總誤差;

(ii)在給定的簇上面進(jìn)行k均值聚類;

(iii)計算將該簇一分為二之后的總誤差;

(3)選擇誤差最小的那個簇進(jìn)行劃分操作.

5 實驗結(jié)果分析與方法驗證

5.1 基于最小二乘法的曲線擬合結(jié)果及分析

基于2017年3月5日到2017年3月17日之間一天24 h的13組統(tǒng)計數(shù)據(jù)(每組數(shù)據(jù)有24個統(tǒng)計點),按一天24個統(tǒng)計點,先計算對應(yīng)各統(tǒng)計點每天平均樣本數(shù),再按這組平均數(shù)據(jù)進(jìn)行基于最小二乘法的多項式擬合獲得的趨勢走向如圖4所示.由圖中連續(xù)的擬合曲線可以看出,一天之中的高峰時段分布在UTC時間的3點、9點和15點內(nèi),圖中離散點為2017年3月5日到2017年3月17日之間,24個統(tǒng)計點所對應(yīng)的實測統(tǒng)計值.

經(jīng)計算對應(yīng)時間點的平均樣本數(shù)的相關(guān)性,以r=0.5為閾值,其中相關(guān)系數(shù)在閾值以上的占92%,因此我們認(rèn)為得到的高峰時段的數(shù)據(jù)是可靠的,基于此,進(jìn)行航跡聚類的實驗.

圖4 2017年3月5日—2017年3月17日之間每日平均樣本數(shù)擬合分布趨勢Fig.4 Fitting of the daily average sample counting distribution between 2017-03-05—2017-03-17

5.2 高峰時段航跡分布范圍聚類劃分

采用實測數(shù)據(jù)對航跡聚類算法進(jìn)行測試,數(shù)據(jù)分為6組,前3組作為測試劃分范圍,后3組驗證方法的有效性.第1組為3月5日–3月11日UTC時間3點內(nèi)所接收的全部航跡點,第2組為3月5日–3月11日UTC時間9點內(nèi)所接收的全部航跡點,第3組為3月5日–3月11日UTC時間15點內(nèi)所接收的全部航跡點,第4組為3月12日–3月17日UTC時間3點內(nèi)所接收的全部航跡點,第5組為3月12日–3月17日UTC時間9點內(nèi)所接收的全部航跡點,第6組為3月12日–3月17日UTC時間15點內(nèi)所接收的全部航跡點,由1、2、3組航跡數(shù)據(jù)聚類結(jié)果分別如圖5(左側(cè))、6(左側(cè))、7(左側(cè))所示.根據(jù)1、2、3組數(shù)據(jù)航跡聚類結(jié)果我們可以劃定高峰時段航跡出現(xiàn)的主要范圍區(qū)域.圖中紅色圓點代表新疆天文臺總部大樓,藍(lán)色圓點代表航跡點,紫色、紅色、綠色圈定的區(qū)域為在此時段內(nèi)航跡點主要聚集范圍.

圖5 UTC 3點航跡聚類(左)和驗證航跡聚類(右)Fig.5 Track clustering during UTC 3(left)and track clustering verifying during UTC 3(right)

5.3 方法驗證

根據(jù)實驗劃定高峰時段范圍區(qū)域,將后3組高峰時段航跡點作為驗證數(shù)據(jù),聚類結(jié)果如圖5(右側(cè))、6(右側(cè))、7(右側(cè))所示.由圖中可以看出,后3組航跡點經(jīng)聚類后多數(shù)落入劃定區(qū)域范圍內(nèi).統(tǒng)計后3組高峰時段航跡點出現(xiàn)在劃定范圍區(qū)域的概率,我們將高峰時段接收的總的航跡點的集合稱作S={s1,s2···sn},落在區(qū)域范圍外的航跡點的集合稱為落在區(qū)域范圍內(nèi)概率公式為:驗證測試所得結(jié)果如表4.經(jīng)過聚類算法的處理結(jié)果表明:航跡在高峰時段的范圍比較確定,符合實際航跡分布特點.

圖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 概率統(tǒng)計Table 4 Probability statistics

6 結(jié)論

以新疆天文臺本部大樓為測試點,針對在測試點接收解析的時間、經(jīng)緯度的航跡樣本數(shù)據(jù),對數(shù)據(jù)按照時間節(jié)點初步統(tǒng)計篩選后,根據(jù)數(shù)據(jù)分布的時間特點,選用基于最小二乘法的多項式對數(shù)據(jù)進(jìn)行擬合,獲得航跡樣本在一天的分布趨勢,經(jīng)計算相關(guān)系數(shù),驗證了用最小二乘法的多項式分析趨勢的可靠性.

根據(jù)高峰時段的分布,將獲得的航跡樣本點分為兩組,一組用于測試,另一組用于驗證方法的有效性.在對第一組高峰時段數(shù)據(jù)進(jìn)行處理之后獲得高峰時段航跡樣本主要分布范圍后,利用相同的方法對第二組數(shù)據(jù)進(jìn)行航跡樣本的聚類處理,經(jīng)驗證航跡樣本點高達(dá)96%落在劃定的區(qū)域范圍內(nèi),驗證了該方法的有效性及實用性.通過數(shù)據(jù)擬合可有效地分析一天之中的高峰時段的分布,根據(jù)擬合結(jié)果對高峰時段航跡進(jìn)行處理,劃定航跡出現(xiàn)的主要范圍區(qū)域,即可獲得一天之后高峰時段的航跡分布范圍,從而為射電天文觀測時段提供有效的參考,最大可能地避開射電天文觀測中CAS對射電天文觀測的影響和對射電望遠(yuǎn)鏡的危害,為規(guī)劃建設(shè)中的射電天文臺址劃分寧靜區(qū)提供與相關(guān)部門交涉的依據(jù).

[1]苗可可,王壯,程翥,等.天文研究與技術(shù),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]張群濤,劉奇,孫正文,等.科學(xué)技術(shù)與工程,2016,1:191

[5]李耀華.基于LMS的射電天文大天線抗干擾技術(shù)的研究與實現(xiàn).長沙:國防科技大學(xué),2011:8

[6]吳忠祖,Haynes M P,Giovanlli R,等.天文學(xué)報,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