曲靖非相干散射雷達(dá)空間碎片微動(dòng)特征提取方法

孫大地 周珂 任同燦

（1.河南師范大學(xué)國際教育學(xué)院, 新鄉(xiāng) 453007；2.河南大學(xué)計(jì)算機(jī)與信息工程學(xué)院, 開封 475004；3.河南省大數(shù)據(jù)分析與處理重點(diǎn)實(shí)驗(yàn)室, 開封 475004）

0 引 言

空間碎片是“地球軌道上或重返大氣層的無功能人造物體，包括其殘塊或部件”[1].目前在軌的空間碎片數(shù)量十分龐大，據(jù)歐空局空間碎片統(tǒng)計(jì)模型估算，超過10 cm 的碎片約有3.6 萬個(gè)，尺寸大小在1 mm至10 cm 之間的碎片更是數(shù)以億計(jì)，對(duì)航天器的在軌運(yùn)行危害巨大，獲取空間碎片運(yùn)動(dòng)特征對(duì)維護(hù)空間活動(dòng)安全和空間的可持續(xù)性利用意義重大[2-5].空間碎片姿態(tài)不受控，在空間攝動(dòng)力(重力梯度、太陽光壓)及殘余角動(dòng)量作用下除具有質(zhì)心平動(dòng)之外，常表現(xiàn)出復(fù)雜的小幅運(yùn)動(dòng)形式(單自旋、平旋或自由翻滾等)，即微動(dòng)特征，它是目前空間目標(biāo)特性研究與目標(biāo)識(shí)別的重要特征量，是建立空間環(huán)境模型的基礎(chǔ)性參數(shù)[6].

非相干散射雷達(dá)是目前地面監(jiān)測電離層的強(qiáng)大手段，由于其具有發(fā)射功率大(一般為MW 量級(jí))、天線增益高(一般大于40 dB)、接收系統(tǒng)噪聲溫度低等技術(shù)特點(diǎn)，在空間碎片探測等方面具有重要應(yīng)用前景.歐洲非相干散射科學(xué)聯(lián)合會(huì)(European Incoherent Scatter Scientific Association, EISCAT)、美國、俄羅斯等利用非相干散射雷達(dá)開展了大量空間碎片探測與研究工作[7-8]，為全球空間碎片監(jiān)測與建模提供了重要支撐.在國家子午工程支持下，中國電波傳播研究所于2012 年在云南曲靖建成我國首套非相干散射雷達(dá)[9].金旺等學(xué)者的研究了表明該雷達(dá)用于空間碎片測量的可行性，系統(tǒng)升級(jí)后有望探測5 cm 左右的空間碎片[10].近年來研究人員利用該雷達(dá)開展了空間碎片探測實(shí)驗(yàn)研究，提取了空間碎片距離、速度、雷達(dá)散射截面(radar cross section, RCS)等主要參數(shù)，獲取了曲靖地區(qū)空間碎片的分布特征[11-14].楊嵩等學(xué)者提出了一種基于窄帶雷達(dá)低頻區(qū)RCS 序列的空間目標(biāo)尺寸和形狀反演方法，并通過理論分析驗(yàn)證了其在曲靖非相干散射雷達(dá)用于空間碎片監(jiān)測的有效性[15].Jian Hu 等學(xué)者提出了一種基于窄帶低分逆合成孔徑雷達(dá)空間碎片的特征提取方法，并仿真驗(yàn)證了方法的有效性和魯棒性[16].這些研究表明地基雷達(dá)是空間碎片探測的有效手段，特別對(duì)曲靖非相干散射雷達(dá)在空間碎片微動(dòng)特征獲取能力進(jìn)行探索意義重大.

隨著人類空間與航天活動(dòng)的不斷增加，在地球軌道上的空間碎片數(shù)目日益增多，加強(qiáng)空間碎片的探測與預(yù)警成為國際上的普遍共識(shí).開展空間碎片微動(dòng)特征提取，有助于認(rèn)識(shí)空間碎片的精細(xì)運(yùn)動(dòng)特性，提升空間碎片監(jiān)測與空間環(huán)境感知能力，同時(shí)進(jìn)一步發(fā)揮曲靖非相干散射雷達(dá)的綜合效益.本文嘗試設(shè)計(jì)一種雷達(dá)目標(biāo)微動(dòng)特征提取流程并應(yīng)用于曲靖非相干散射雷達(dá)，同時(shí)通過仿真和實(shí)測數(shù)據(jù)對(duì)所提流程進(jìn)行了分析驗(yàn)證.

1 微動(dòng)特征提取方法

1.1 微動(dòng)目標(biāo)運(yùn)動(dòng)模型

Chen V.C.首次將微多普勒的概念引入到雷達(dá)目標(biāo)識(shí)別領(lǐng)域，并推導(dǎo)了振動(dòng)、轉(zhuǎn)動(dòng)、錐動(dòng)以及翻滾四種基本的微動(dòng)數(shù)學(xué)模型[17]，指出利用目標(biāo)的微多普勒頻率可反演目標(biāo)的運(yùn)動(dòng)特征(如轉(zhuǎn)動(dòng)半徑和轉(zhuǎn)動(dòng)頻率等)，通過時(shí)頻變換進(jìn)行雷達(dá)成像可反演目標(biāo)的結(jié)構(gòu)參數(shù)(即散射中心在水平面上的投影分布).

在RSW 坐標(biāo)系下微動(dòng)目標(biāo)運(yùn)動(dòng)如圖1 所示.由于RSW 坐標(biāo)系的三軸分別為目標(biāo)的徑向方向、沿跡方向以及垂跡方向，故目標(biāo)微動(dòng)中心在該坐標(biāo)系下相對(duì)雷達(dá)的方位角和俯仰角均可視為0.假設(shè)目標(biāo)自轉(zhuǎn)軸與垂跡方向重合，則微動(dòng)目標(biāo)的運(yùn)動(dòng)模型可簡化為圖1 下圖所示.目標(biāo)回波相位表達(dá)式φm-D(t)為

圖1 微動(dòng)目標(biāo)在RSW 坐標(biāo)系下的運(yùn)動(dòng)示意圖(上)及其簡化情況(下)Fig.1 The motion diagram(top) and its simplifying(bottom)of the microtremor target in the RSW coordinate

式中：R0為 微動(dòng)中心距雷達(dá)的初始距離；rW和rS分別為目標(biāo)微動(dòng)半徑(自旋/翻滾)在垂直面(R-W)和水平面(R-S)上的投影； ωW和 ωS分別為目標(biāo)的轉(zhuǎn)動(dòng)軸在垂跡和沿跡方向上的分量； θ0和 ?0分別為散射中心和轉(zhuǎn)動(dòng)中心連線與沿跡方向、垂跡方向的初始夾角； λ為雷達(dá)載波波長.

進(jìn)一步地，可得到目標(biāo)的微動(dòng)頻率表達(dá)式：

1.2 微動(dòng)特征提取流程與仿真

雷達(dá)回波經(jīng)過脈壓處理后，目標(biāo)的平動(dòng)運(yùn)動(dòng)得到了補(bǔ)償，但是相位中仍保留有一定的平動(dòng)余項(xiàng)必須修正，否則后續(xù)處理中會(huì)存在多普勒模糊，影響成像質(zhì)量.從回波時(shí)頻變換域的雷達(dá)成像可以直觀提取目標(biāo)的二維幾何形態(tài)和運(yùn)動(dòng)特征，具體提取流程如下：

首先，對(duì)脈壓后的雷達(dá)回波進(jìn)行預(yù)處理，用Keystone 變換或搜索出的平動(dòng)運(yùn)動(dòng)參數(shù)修正因目標(biāo)的徑向運(yùn)動(dòng)引入的距離徙動(dòng)，同時(shí)修正平動(dòng)對(duì)雷達(dá)回波產(chǎn)生的相位信息[18].

其次，使用De-chirp 去調(diào)制方法，計(jì)算目標(biāo)平動(dòng)余項(xiàng)的速度和加速度參數(shù)，進(jìn)一步補(bǔ)償回波中余留的平動(dòng)相位調(diào)制，經(jīng)此處理后可認(rèn)為目標(biāo)回波相位只受微動(dòng)的調(diào)制.De-chirp 方法是將輸入信號(hào)與參考信號(hào)相乘運(yùn)算，接著對(duì)產(chǎn)生的信號(hào)進(jìn)行頻譜分析，對(duì)參考信號(hào)的要求是其應(yīng)具有與輸入信號(hào)相同的調(diào)頻率[19].

接著，使用S-method(SM)方法計(jì)算目標(biāo)回波的時(shí)頻分布.SM 方法具有良好的時(shí)頻分辨率且可通過在頻域加窗有效抑制相干干擾項(xiàng)，缺點(diǎn)是對(duì)信噪比要求較高，這是因?yàn)樵摲椒ǖ谋举|(zhì)是利用目標(biāo)回波能量分布來計(jì)算其時(shí)頻分布[20-22].當(dāng)目標(biāo)信噪比較低時(shí)，可對(duì)回波進(jìn)行短時(shí)間傅里葉變換，利用時(shí)頻譜的幅度和相位的相干性實(shí)現(xiàn)對(duì)目標(biāo)時(shí)頻分布的高分辨率聚焦，本文僅考慮高信噪比條件下的目標(biāo)微動(dòng)特征提取.

最后，通過Inverse-Radon 變換(inverse Radon transform, IRT)將時(shí)頻像投影到二維空間[23].IRT 具有將正弦曲線映射為參數(shù)空間中一點(diǎn)的變換特性，因此可以對(duì)微動(dòng)目標(biāo)進(jìn)行頻域成像.由式(2)可知微動(dòng)目標(biāo)的正弦調(diào)頻(sinusoid frequency modulation,SFM)信號(hào)是關(guān)于時(shí)間的正弦函數(shù)，因此IRT 可以將單個(gè)散射中心的時(shí)頻能量聚焦到變換域上的一個(gè)點(diǎn)，具有很好的魯棒性.

設(shè)散射中心i的SFM 信號(hào)為

需要指出的是，因?yàn)檗D(zhuǎn)動(dòng)頻率 ω為未知量，因此IRT 的弧度自變量也為未知量，故需要對(duì)可能的轉(zhuǎn)動(dòng)頻率逐一代入計(jì)算.由于上一步De-chirp 處理時(shí)已遍歷搜索得到目標(biāo)的微動(dòng)頻率初值，此時(shí)僅需在該頻率及其倍數(shù)附近搜索即可(即考慮微多普勒模糊).當(dāng)搜索的頻率與真值相接近時(shí)，IRT 后的能量將十分集中(聚焦到一點(diǎn))，通過下式可計(jì)算IRT 的集中度(MC).

式 中，R(x,y)為 IRT 的 結(jié) 果.將 具 有 最 大MC值 的IRT 對(duì)應(yīng)周期數(shù)作為后驗(yàn)值代入，并除以 2ω/λ，可以得到最終真實(shí)距離坐標(biāo)下的散射中心在二維空間上的分布.

圖2 給出了空間微動(dòng)目標(biāo)4 個(gè)散射中心的仿真結(jié)果.在仿真中使用Born 一階近似理論描繪目標(biāo)的散射回波過程，假設(shè)各散射中心的轉(zhuǎn)動(dòng)半徑分別為8、3、3 和8 個(gè)波長，初始相角分別為45°、135°、-45°與-135°，且該目標(biāo)具有3 m/s 的平動(dòng)速度余項(xiàng)與0.5 m/s2的平動(dòng)加速度余項(xiàng)，將以上參數(shù)代入式(1)求出仿真目標(biāo)的回波相位變化，進(jìn)而可仿真出目標(biāo)的回波信號(hào).圖2(a)為回波信號(hào)的原始時(shí)頻分布，可以看出SFM 信號(hào)偏離零頻軸，無法直接進(jìn)行IRT.圖2(b)為對(duì)回波信號(hào)進(jìn)行De-chirp 的處理過程，其中第一行為第1 時(shí)延(lag-1)共軛時(shí)延積的快速傅里葉變換(fast Fourier transform, FFT)結(jié)果，右側(cè)為FFT 最大響應(yīng)部分局部放大，紅線代表真值(2aτ/λ)，可以看出FFT 最大響應(yīng)頻率近似等于真值；第二行為第n時(shí)延(lag-n)共軛時(shí)延積的FFT 結(jié)果，從右側(cè)的局部放大可以看出FFT 最大響應(yīng)頻率近似等于真值(4v/λ)；第三行代表lag-n時(shí)延積的FFT 的熵隨時(shí)延的變化，可以看出當(dāng)時(shí)延等于半個(gè)轉(zhuǎn)動(dòng)周期時(shí)，F(xiàn)FT 具有最小熵.圖2(c)為補(bǔ)償目標(biāo)平動(dòng)余項(xiàng)后的時(shí)頻圖，可以看出此時(shí)各散射中心的SFM 信號(hào)(黑色虛線)均關(guān)于零頻軸對(duì)稱，且只顯現(xiàn)出正弦調(diào)制.圖2(d)為對(duì)圖2(c)時(shí)頻圖進(jìn)行IRT 后的結(jié)果，可以看出4 個(gè)散射中心在IRT 后分別聚焦到不同的位置，對(duì)應(yīng)于散射中心的起始位置.白色“+”代表散射中心的真實(shí)位置，可以看出成像位置與真實(shí)位置基本重合，經(jīng)計(jì)算得到各散射中心的平均測量偏差為0.2λ(偏差主要來自計(jì)算時(shí)頻分布時(shí)的交叉項(xiàng)影響，圖中近似波紋的形狀)，驗(yàn)證了IRT 用于提取目標(biāo)微動(dòng)頻率以及各散射中心位置(轉(zhuǎn)動(dòng)半徑)的有效性.

圖2 四散射中心空間目標(biāo)仿真結(jié)果Fig.2 The simulation of the space microtremor target with 4 scattering centers

2 實(shí)測數(shù)據(jù)分析

利用曲靖非相干散射雷達(dá)開展實(shí)驗(yàn)驗(yàn)證，實(shí)驗(yàn)期間雷達(dá)發(fā)射波形為13 位巴克碼，脈沖重復(fù)周期為12 ms，碼元寬度為30 μs，回波采樣帶寬為6.25 MHz，探測模式為凝視探測，利用上文方法對(duì)采集的原始回波數(shù)據(jù)進(jìn)行處理，處理過程中暫未考慮碎片在波束內(nèi)劃過時(shí)在波束中不同位置強(qiáng)度衰減的影響.

圖3 為2017-06-30T04:55:33UT 某空間碎片(軌道高度約448.5 km)的分析結(jié)果.圖3(a)為該空間碎片目標(biāo)回波的原始時(shí)頻圖，可見頻譜能量集中分布在一SFM 信號(hào)附近，且該曲線不顯現(xiàn)明顯的正弦調(diào)制.這是由于在脈壓時(shí)沒有完全補(bǔ)償平動(dòng)運(yùn)動(dòng)，致使SFM 信號(hào)呈現(xiàn)出直線疊加正、余弦曲線的形態(tài).利用De-chirp 方法找到目標(biāo)的平動(dòng)速度余項(xiàng)和加速度余項(xiàng)后，對(duì)平動(dòng)運(yùn)動(dòng)引入的相位調(diào)制進(jìn)行補(bǔ)償，再利用SM 方法重新計(jì)算目標(biāo)回波的時(shí)頻分布，得到圖3(b)結(jié)果.為便于下一步IRT，將時(shí)頻圖在多普勒維上進(jìn)行了歸一化處理(利用各時(shí)間上頻譜的最大值).在IRT 時(shí)需要遍歷搜索微動(dòng)周期個(gè)數(shù)，通過測量相應(yīng)IRT 結(jié)果的集中度，找到MC值最高的IRT 作為最終的變換結(jié)果，此時(shí)對(duì)應(yīng)的周期數(shù)可視為時(shí)頻分布中含有完整微動(dòng)周期的個(gè)數(shù).不同周期數(shù)對(duì)應(yīng)的MC結(jié)果如圖3(c)所示，可以看出最大MC值對(duì)應(yīng)1 個(gè)周期，與圖3(b)中SFM 信號(hào)周期的個(gè)數(shù)相一致.最終的IRT 成像結(jié)果如圖3(e)所示，可以看出該碎片目標(biāo)具有一個(gè)強(qiáng)散射中心(周圍的兩個(gè)散射中心可能來自能量泄露)和一個(gè)弱散射中心.將各散射中心及其轉(zhuǎn)動(dòng)中心在極坐標(biāo)系下進(jìn)行排列，得到圖3(d)中的結(jié)果，可以看出該目標(biāo)形狀大體為長條形、兩散射中心的初始姿態(tài)角分別為-55°和-76°、翻滾半徑分別約為2.24 m 和2.88 m、翻滾頻率和周期分別為0.08 Hz和12 s.通過比較時(shí)頻能量的變化，可以看出目標(biāo)在2 s 時(shí)刻具有最強(qiáng)回波功率，此時(shí)散射中心姿態(tài)角約為0°；而在8～9 s 時(shí)回波功率驟降，此時(shí)姿態(tài)角約為180°，這反映了該目標(biāo)散射中心姿態(tài)的變化，也間接反映了該目標(biāo)形狀的不規(guī)則性.

圖3 2017-06-30T04:55:33UT 某空間碎片實(shí)驗(yàn)分析結(jié)果Fig.3 The experimental result of space debris at 04:55:33UT on June 30, 2017

圖4 為2015-01-17T03:49:05UT 某空間碎片(軌道高度約1 012.8 km)的分析結(jié)果.可以看出該目標(biāo)時(shí)頻分布較復(fù)雜，具有多個(gè)散射中心，IRT 成像結(jié)果也印證了這一點(diǎn)，從圖4(e)可以看出該目標(biāo)具有5 個(gè)散射中心.其中最強(qiáng)散射中心的微動(dòng)參數(shù)為：翻滾半徑0.4 m、翻滾頻率0.35 Hz、翻滾周期2.83 s.該目標(biāo)的時(shí)頻圖中可能存在干擾項(xiàng)，而且各周期數(shù)自變量對(duì)應(yīng)的MC值較接近，尚需進(jìn)一步分析.

圖4 2015-01-17T03:49:05UT 某空間碎片實(shí)驗(yàn)分析結(jié)果(子圖含義與圖3 一致)Fig.4 The experimental results of space debris on January 17, 2015 (the content of different sub-graphs is the same as Fig.3 )

圖5 為2017-06-26T11:41:27UT 某空間碎片(軌道高度約1 180.2 km)的分析結(jié)果.可看出其自旋半徑約6.5 m、自旋頻率約0.05 Hz、自旋周期約18.9 s，初步認(rèn)為其具有自旋穩(wěn)定的運(yùn)動(dòng)姿態(tài).該目標(biāo)尺寸較大且姿態(tài)角在過境期間變化較大(約228°)，IRT 圖像顯示該目標(biāo)具有十分接近的散射中心，初步分析是因?yàn)槠湫螤畈灰?guī)則所致.

圖5 2017-06-26T11:41:27UT 某空間碎片實(shí)驗(yàn)分析結(jié)果(子圖含義與圖3 一致)Fig.5 The experimental results of space debris at 11:41:27UT on June 26, 2017(the content of different sub-graphs is the same as Fig.3 )

圖6 為2017-06-30T04:16:39UT 某空間碎片(軌道高度約618.6 km)的分析結(jié)果.從時(shí)頻圖分布可猜測該目標(biāo)具有三軸穩(wěn)定的運(yùn)動(dòng)姿態(tài).在補(bǔ)償平動(dòng)余項(xiàng)后，回波信號(hào)的SFM 信號(hào)隨時(shí)間變化十分緩慢.由于穿過雷達(dá)波束期間該目標(biāo)相對(duì)非相干散射雷達(dá)的姿態(tài)角變化很小，故在IRT 像中顯現(xiàn)出單一散射中心，且轉(zhuǎn)動(dòng)半徑很小.由于該目標(biāo)的時(shí)頻分布變化很小，因此其IRT 成像結(jié)果實(shí)際意義不大.

圖6 2017-06-30T04:16:39UT 某空間碎片實(shí)驗(yàn)分析結(jié)果(子圖含義與圖3 一致)Fig.6 The experimental results of space debris at 04:16:39UT on June 30, 2017(the content of different sub-graphs is the same as Fig.3 )

3 結(jié) 論

非相干散射雷達(dá)是目前最強(qiáng)大的電離層地基監(jiān)測手段，在空間碎片監(jiān)測方面具有很大潛力.國內(nèi)外利用非相干散射雷達(dá)開展了大量探測實(shí)驗(yàn)工作，提取了空間碎片距離、速度、RCS 等主要參數(shù)，獲取了不同地區(qū)的空間碎片分布特征，為空間碎片模型研究和預(yù)警服務(wù)提供了重要數(shù)據(jù)支撐，但是對(duì)空間碎片目標(biāo)特征提取工作，包括微動(dòng)特征、幾何與結(jié)構(gòu)參數(shù)等的提取尚未深入開展.

本文借鑒雷達(dá)目標(biāo)微動(dòng)特征提取方法，嘗試從非相干散射雷達(dá)實(shí)測數(shù)據(jù)中提取了空間碎片微動(dòng)特征，包括空間碎片翻滾半徑、周期、初始姿態(tài)角、衛(wèi)星自旋半徑與周期等.實(shí)測數(shù)據(jù)與理論仿真基本一致，且符合實(shí)際特征規(guī)律，驗(yàn)證了利用非相干散射雷達(dá)開展空間碎片微動(dòng)特征提取的可行性.

本文方法認(rèn)為目標(biāo)轉(zhuǎn)動(dòng)軸與垂跡方向重合，IRT 成像是將目標(biāo)各散射點(diǎn)投影到水平面上的結(jié)果，下一步需考慮這種簡化對(duì)目標(biāo)特征提取的影響(即修正由轉(zhuǎn)動(dòng)軸-視線夾角引入的縮放效應(yīng)).此外需利用更多數(shù)據(jù)，分析該方法在不同情況下的有效性(如信噪比條件、多散射中心、時(shí)頻圖中存在交叉干擾項(xiàng)等)，并與其他方法相對(duì)比(如SRDI 單距離門-多普勒干涉成像)，深入探討高效的空間碎片微動(dòng)特征提取方法.

致謝：感謝中國電子科技集團(tuán)公司第二十二研究所、云南昆明電磁波環(huán)境國家野外科學(xué)觀測研究站提供曲靖非相干散射雷達(dá)實(shí)驗(yàn)數(shù)據(jù).