陳是扦，彭志科，邢冠培

(1.上海交通大學(xué) 機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200240； 2. 上海航天電子技術(shù)研究所，上海 201109)

0 引言

雷達(dá)通過發(fā)射、接收目標(biāo)反射的電磁波進(jìn)行探測(cè)。當(dāng)目標(biāo)以一定速度運(yùn)動(dòng)時(shí)，其反射信號(hào)將偏離雷達(dá)發(fā)射信號(hào)的頻率，該現(xiàn)象稱為多普勒效應(yīng)。如果目標(biāo)除了整體平移運(yùn)動(dòng)外，還進(jìn)行局部振動(dòng)、轉(zhuǎn)動(dòng)等微運(yùn)動(dòng)，則會(huì)使雷達(dá)回波信號(hào)產(chǎn)生頻率調(diào)制，并在目標(biāo)整體平移運(yùn)動(dòng)引起的多普勒頻率周圍產(chǎn)生邊帶，這種現(xiàn)象稱為微多普勒效應(yīng)[1-3]。微多普勒調(diào)制頻率可反映目標(biāo)的微動(dòng)特性和結(jié)構(gòu)特征，被廣泛應(yīng)用于雷達(dá)目標(biāo)分類識(shí)別中[4-5]。旋轉(zhuǎn)運(yùn)動(dòng)是航天領(lǐng)域中最為常見的微運(yùn)動(dòng)，如衛(wèi)星天線轉(zhuǎn)動(dòng)、彈道中段彈頭自旋運(yùn)動(dòng)等，因此本文將主要關(guān)注旋轉(zhuǎn)目標(biāo)的微多普勒頻率。

目標(biāo)運(yùn)動(dòng)使微多普勒信號(hào)頻率具有時(shí)變特征[6]。目前，國(guó)內(nèi)外已有大量利用時(shí)頻分析方法處理這類非平穩(wěn)信號(hào)的報(bào)道,如：文獻(xiàn)[7]利用經(jīng)驗(yàn)?zāi)Ｊ椒纸夥椒ǚ蛛x出目標(biāo)旋轉(zhuǎn)部件產(chǎn)生的微多普勒特征分量，以消除微多普勒特征對(duì)目標(biāo)主體部分ISAR(逆合成孔徑雷達(dá))成像的影響，但經(jīng)驗(yàn)?zāi)Ｊ椒纸夥椒ㄈ菀桩a(chǎn)生邊界效應(yīng)、模式混疊等問題，難以處理航天領(lǐng)域各種復(fù)雜的雷達(dá)信號(hào)；文獻(xiàn)[8-9]分別利用Fourier-Bessel級(jí)數(shù)和Chirplet變換提取雷達(dá)回波信號(hào)中的微多普勒特征，利用Fourier-Bessel級(jí)數(shù)分離雷達(dá)信號(hào)中的非平穩(wěn)成分與平穩(wěn)成分，而非平穩(wěn)成分往往與雷達(dá)微多普勒特征相對(duì)應(yīng)，不能進(jìn)一步分離不同散射點(diǎn)的微多普勒特征。上述文獻(xiàn)主要介紹了微多普勒信號(hào)分解，未見信號(hào)頻率特征提取的報(bào)道。文獻(xiàn)[10-11]通過檢測(cè)回波信號(hào)時(shí)頻分布的局部峰值點(diǎn)，估計(jì)目標(biāo)的微多普勒頻率，但由于目標(biāo)不同散射點(diǎn)的微多普勒頻率分量通常在時(shí)頻分布上相互重疊，該方法難以得到感興趣的頻率分量，且信號(hào)時(shí)頻分布的分辨率有限，因此較難估計(jì)高精度的微多普勒頻率；文獻(xiàn)[12]提出一種參數(shù)化解調(diào)方法，用于估計(jì)多項(xiàng)式頻率調(diào)制信號(hào)的瞬時(shí)頻率，而該方法能有效處理瞬時(shí)頻率交叉、重疊的信號(hào)，但其采用的多項(xiàng)式模型在逼近復(fù)雜瞬時(shí)頻率曲線時(shí)存在龍格現(xiàn)象，無法估計(jì)高度振蕩的微多普勒頻率。

鑒于旋轉(zhuǎn)目標(biāo)微多普勒信號(hào)具有正弦頻率調(diào)制特征，本文提出了一種基于正弦模型的參數(shù)化解調(diào)方法。該方法與文獻(xiàn)[12]中的多項(xiàng)式模型解調(diào)方法不同，將頻率估計(jì)問題轉(zhuǎn)化為模型參數(shù)優(yōu)化問題，利用粒子群算法優(yōu)化參數(shù)，對(duì)旋轉(zhuǎn)目標(biāo)呈周期振蕩的微多普勒頻率進(jìn)行有效估計(jì)。

1 理論背景

1.1 信號(hào)模型

將旋轉(zhuǎn)目標(biāo)簡(jiǎn)化為繞固定旋轉(zhuǎn)中心運(yùn)動(dòng)的散射點(diǎn)[13]，雷達(dá)與旋轉(zhuǎn)散射點(diǎn)的幾何關(guān)系如圖1所示。圖中：散射點(diǎn)P繞回轉(zhuǎn)中心O以角頻率fr逆時(shí)針旋轉(zhuǎn)；l0為雷達(dá)到回轉(zhuǎn)中心的距離；l1為回轉(zhuǎn)半徑；LP為雷達(dá)到散射點(diǎn)的瞬時(shí)距離；φ0、φt分別為起始時(shí)刻和t時(shí)刻散射點(diǎn)P的旋轉(zhuǎn)角，φt=φ0+2πfrt。

圖1 雷達(dá)與旋轉(zhuǎn)散射點(diǎn)的幾何關(guān)系Fig.1 Geometrical relationship between radar and rotating scattering point

由圖1中的幾何關(guān)系得到雷達(dá)與散射點(diǎn)之間距離，即

(1)

假設(shè)散射體位于雷達(dá)遠(yuǎn)場(chǎng)，即(l1/l0)2→0，則式(1)可簡(jiǎn)化為

LP(t)≈l0+l1cos(φ0+2πfrt)

(2)

因此，點(diǎn)P的雷達(dá)回波基帶信號(hào)可表示為

(3)

(4)

(5)

式中：K為散射點(diǎn)數(shù)目；σi為第i個(gè)散射點(diǎn)的散射系數(shù)；Li(t)為雷達(dá)與第i個(gè)散射點(diǎn)之間的瞬時(shí)距離。

假設(shè)所有微動(dòng)散射點(diǎn)的旋轉(zhuǎn)頻率為fr，則參照式(2)可得到瞬時(shí)距離函數(shù)，其表達(dá)式為

Li(t)≈l0+li1cos(φi0+2πfrt)

(6)

式中：li1為第i個(gè)散射點(diǎn)的旋轉(zhuǎn)半徑；φi0為第i個(gè)散射點(diǎn)的起始旋轉(zhuǎn)角。參照式(4)可得到第i個(gè)散射點(diǎn)的微多普勒頻率，即

(7)

由式(7)可得：不同散射點(diǎn)的微多普勒頻率是一些幅值和相位不同、角頻率相同的正弦函數(shù)，這些微多普勒頻率通常在時(shí)頻面上相互重疊，難以提取。

1.2 參數(shù)化解調(diào)

參數(shù)化解調(diào)方法可估計(jì)特定信號(hào)的模型參數(shù)。廣義調(diào)頻信號(hào)模型為

(8)

式中：A為信號(hào)幅值；ψ0為信號(hào)初始相位；fI(t)為信號(hào)瞬時(shí)頻率。若式(8)為雷達(dá)回波信號(hào)，則fI(t)為微多普勒頻率，其表達(dá)式見式(4)。

參數(shù)化解調(diào)算子的表達(dá)式為

Φ(t；P,fc)=

(9)

式中：fc為載波頻率；g(t;P)為特定的參數(shù)化解調(diào)函數(shù)[12]；P為解調(diào)函數(shù)的參數(shù)集合。將式(8)乘以式(9)得到解調(diào)信號(hào)，即

(10)

當(dāng)且僅當(dāng)fI(t)=g(t;P)時(shí)，式(10)為單頻信號(hào)sd(t)=Aexp[j(2πfct+ψ0)]，即信號(hào)能量集中在頻率fc處，信號(hào)頻譜值在fc處達(dá)到最大。因此，可通過優(yōu)化信號(hào)模型參數(shù)P使sd(t)在載波頻率fc處的頻譜值達(dá)到最大，即

(11)

需要注意的是，原始信號(hào)由解調(diào)信號(hào)調(diào)制得到，即

s(t)=sd(t)×conj[Φ(t;P,fc)]=

(12)

式中：conj為復(fù)數(shù)共軛。

2 基于參數(shù)化解調(diào)的微多普勒頻率估計(jì)方法

由于旋轉(zhuǎn)目標(biāo)的微多普勒頻率隨時(shí)間正弦變化，因此提出基于正弦模型的參數(shù)化解調(diào)方法，用于提取目標(biāo)微多普勒頻率。將式(7)整理為

(13)

式中：ai為多普勒頻率幅值，ai=4πfrli1/λ；φi0=φi0+π。參照式(9)，構(gòu)造與式(13)相匹配的參數(shù)化解調(diào)算子，即

(14)

式中：fc為載波頻率，fc=0(參考式(9))；Pi={ai,fr,φi0}為待估計(jì)的模型參數(shù)集。參照式(11)可估計(jì)微多普勒頻率參數(shù)，即

(15)

式中：sr(t)為雷達(dá)回波信號(hào)(見式(5))。采用粒子群優(yōu)化算法對(duì)式(15)進(jìn)行優(yōu)化求解[14]。

回波信號(hào)包含K個(gè)信號(hào)分量，由于強(qiáng)信號(hào)分量對(duì)弱信號(hào)分量的估計(jì)有抑制作用，通常只能估計(jì)到能量最強(qiáng)(即σi最大)的信號(hào)分量參數(shù)。為得到所有K個(gè)散射點(diǎn)的微多普勒頻率參數(shù)，采用迭代估計(jì)算法，將估計(jì)到的強(qiáng)信號(hào)分量從原始信號(hào)中剔除(消除強(qiáng)信號(hào)分量的干擾)，從剩余信號(hào)中估計(jì)其他的弱信號(hào)分量參數(shù)。

算法流程具體步驟如下：

1) 令i=1；

7) 重復(fù)執(zhí)行步驟2)～6)，直到i≥K。

上述步驟中，步驟6)將當(dāng)前提取到的強(qiáng)信號(hào)分量從回波信號(hào)中移除，因此該算法能在下次迭代中從剩余信號(hào)估計(jì)到其他能量較弱的信號(hào)分量參數(shù)。上述算法經(jīng)過K次迭代可提取所有散射點(diǎn)的微多普勒頻率參數(shù)。算法流程如圖2所示。

圖2 算法流程Fig.2 Flowchart of algorithm

3 算法驗(yàn)證

3.1 仿真驗(yàn)證

雷達(dá)波長(zhǎng)λ=0.06 m；雷達(dá)與旋轉(zhuǎn)目標(biāo)中心的距離l0=500 m；目標(biāo)包含3個(gè)強(qiáng)散射點(diǎn)，其旋轉(zhuǎn)半徑分別為l11=3 m，l21=2.5 m，l31=1 m；初始旋轉(zhuǎn)角分別為φ10=π/3 rad，φ20=5π/6 rad，φ30=4π/3 rad；散射系數(shù)σ1=3，σ2=2.8，σ3=2；目標(biāo)旋轉(zhuǎn)角頻率fr=2 Hz。假設(shè)信號(hào)被高斯白噪聲污染，信噪比為0 dB，則利用短時(shí)傅里葉變換得到含噪聲的雷達(dá)回波信號(hào)的時(shí)頻分布，仿真信號(hào)時(shí)頻分布如圖3所示。圖中可見，3個(gè)散射點(diǎn)的微多普勒頻率曲線在時(shí)頻面上相互相交、重疊。

圖3 仿真信號(hào)時(shí)頻分布Fig.3 Time-frequency distribution of simulated signal

采用本文方法估計(jì)式(13)中的微多普勒模型參數(shù)，包括微多普勒振幅ai、振蕩頻率fr(目標(biāo)旋轉(zhuǎn)角頻率)、初始相位φi0。進(jìn)一步估計(jì)散射點(diǎn)的旋轉(zhuǎn)半徑li1，li1=λai/(4πfr)。3個(gè)散射點(diǎn)的參數(shù)估計(jì)結(jié)果見表1。結(jié)果表明：本文方法在強(qiáng)噪聲干擾下能得到高精度的參數(shù)估計(jì)結(jié)果。

表1 仿真信號(hào)微多普勒參數(shù)估計(jì)結(jié)果

圖4 仿真信號(hào)微多普勒頻率提取結(jié)果對(duì)比Fig.4 Comparison of extracted micro-Doppler frequencies of simulated signal

分別用文獻(xiàn)[11]中的方法和本文方法進(jìn)行微多普勒頻率提取，對(duì)比結(jié)果如圖4所示。文獻(xiàn)[11]通過檢測(cè)回波信號(hào)時(shí)頻分布(見圖3)的局部峰值點(diǎn)提取微多普勒頻率曲線。由圖4可見：不同曲線在交點(diǎn)處模糊，局部峰值檢測(cè)法在通過微多普勒頻率交叉點(diǎn)后容易識(shí)別到錯(cuò)誤信號(hào)分量，而本文方法能正確識(shí)別交叉信號(hào)分量。

3.2 實(shí)驗(yàn)驗(yàn)證

為驗(yàn)證本文方法在實(shí)際應(yīng)用中的有效性，對(duì)雷達(dá)實(shí)測(cè)數(shù)據(jù)的微多普勒頻率進(jìn)行估計(jì)。實(shí)驗(yàn)裝置與文獻(xiàn)[15]中的類似，如圖5所示。圖中：2個(gè)金屬圓筒以180°的間隔布置在轉(zhuǎn)盤上，轉(zhuǎn)盤以一定速度旋轉(zhuǎn)運(yùn)動(dòng)；將2個(gè)金屬圓筒視為2個(gè)旋轉(zhuǎn)目標(biāo)，并用發(fā)射頻率為10.5 GHz(即波長(zhǎng)為0.028 5 m)的連續(xù)波雷達(dá)探測(cè)2個(gè)目標(biāo)的旋轉(zhuǎn)運(yùn)動(dòng)；2個(gè)圓筒距離旋轉(zhuǎn)中心的距離都為0.1 m。

實(shí)測(cè)信號(hào)的時(shí)頻分布如圖6所示。由圖可見：2個(gè)旋轉(zhuǎn)目標(biāo)的微多普勒頻率曲線隨時(shí)間正弦變化，兩曲線相位差約為180°。

用本文方法估計(jì)實(shí)驗(yàn)信號(hào)的微多普勒頻率參數(shù)，結(jié)果見表2。實(shí)驗(yàn)中轉(zhuǎn)盤的實(shí)際轉(zhuǎn)速未知，無法提供各頻率參數(shù)的真實(shí)值，但由估計(jì)到的頻率參數(shù)計(jì)算得到的目標(biāo)旋轉(zhuǎn)半徑li1與實(shí)際半徑(0.1 m)的誤差不超過2%。實(shí)驗(yàn)結(jié)果表明：本文方法對(duì)于實(shí)測(cè)雷達(dá)數(shù)據(jù)仍具有較高的頻率估計(jì)精度。

本文方法與傳統(tǒng)時(shí)頻峰值檢測(cè)法得到的微多普勒頻率曲線對(duì)比如圖7所示。由圖可見：傳統(tǒng)時(shí)頻峰值檢測(cè)法無法正確識(shí)別交叉的信號(hào)分量，且由于其受其他信號(hào)分量的干擾，提取到的頻率曲線出現(xiàn)了局部波動(dòng)和毛刺；本文方法采用與信號(hào)匹配的參數(shù)化模型，抗干擾能力強(qiáng)，能得到精確的微多普勒頻率曲線。

4 結(jié)束語

本文對(duì)旋轉(zhuǎn)目標(biāo)的微多普勒頻率提取問題進(jìn)行了研究，分析了基于等效散射點(diǎn)模型的旋轉(zhuǎn)目標(biāo)雷達(dá)回波信號(hào)解析形式。針對(duì)信號(hào)模型特征，提出了基于正弦模型的參數(shù)化解調(diào)方法，以提取旋轉(zhuǎn)目標(biāo)的微多普勒頻率。在算法具體實(shí)施中，采用迭代頻率提取方法，逐個(gè)獲取每個(gè)散射點(diǎn)的微多普勒頻率。在仿真實(shí)驗(yàn)中，將本文方法與傳統(tǒng)時(shí)頻峰值檢測(cè)法進(jìn)行對(duì)比。結(jié)果表明：本文方法能有效提取相互交叉的微多普勒頻率曲線，算法精度高，抗干擾能力強(qiáng)，得到的微多普勒頻率可進(jìn)一步用于目標(biāo)運(yùn)動(dòng)參數(shù)識(shí)別提取，在識(shí)別、監(jiān)測(cè)彈道導(dǎo)彈飛行參數(shù)等方面具有重要應(yīng)用前景。后續(xù)將基于該方法開展航天領(lǐng)域雷達(dá)目標(biāo)識(shí)別應(yīng)用的系統(tǒng)研究。