基于相位補償的外輻射源雷達旋翼目標特征提取

周子鉑 張朝偉 夏賽強 許道明 高 妍 曾曉雙

①(空軍預警學院 武漢 430019)

②(中國人民解放軍31001部隊 北京 100095)

1 引言

鑒于全球定位導航系統(Global Navigation Satellite System,GNSS)廣泛的覆蓋率和穩定的工作性能，利用GNSS信號進行目標探測與識別的應用已經引起廣泛研究。該型外輻射源雷達無發射機，信號源自空間軌道上的定位導航衛星，而接收機位于地面，屬收發分置雷達系統[1,2]。在目標散射方面具備雷達散射截面增強的特性，即前向散射增強，所以前向散射系數遠高于后向散射系數[3]。此散射特性促使外輻射源雷達在處理GNSS輻射的低功率信號時具備相當的可行性，尤其在近空域目標識別應用之中[4,5]。

不同于目標整體的平動，旋翼葉片的旋轉運動被稱為“微動”[6]。目標的微動導致雷達接收回波中除了含有平動多普勒頻率之外，還包含有“微多普勒”頻率[7,8]。而微多普勒頻率的大小和變化趨勢與旋翼目標的特征密切相關，不同的旋翼目標具有自身獨特的微多普勒效應[9—11]，因此可以利用雷達回波的微多普勒效應提取目標旋翼特征參數。

文獻[12]推導了基于GNSS外輻射源的直升機旋翼目標回波理論模型，詳細分析了微多普勒特性分析對于外輻射源的系統要求，并結合幾種典型直升機的參數特征給出了不同觀測視角下所能達到的作用距離，以及回波信號的信干噪比。此外，Carmine.C等人還利用時頻工具獲取了典型直升機模型的時頻分布。基于單基地模型，朱名爍和陳永彬在文獻[13,14]中分析了葉片上散射點的不同分布情況和目標不同的姿態角對回波的影響，并在此基礎上研究了時域閃爍現象的物理形成機理，但是未進一步提出具體的旋翼目標參數提取方法。針對微多普勒效應分析問題，文獻[15]根據轉動部件的周期性旋轉導致微多普勒頻率呈正弦變化的特點，提出利用霍夫變換檢測時頻圖譜中的微動曲線，進而實現目標參數特征的提取。文獻[16]提出利用正交匹配追蹤(Orthogonal Match Pursuit,OMP)算法在電視廣播信號做輻射源的條件下對旋翼直升機目標進行參數估計，但是匹配追蹤這類貪婪算法搜索到的可能是局域內的極值，而非全域的最值，因此會降低算法的穩定性。張群等在文獻[17]中研究了頻分復用模式下基于多重觀測矢量模型和系數重構方法，該方法可以實現較高的參數估計精度，但是所要求的數據量將隨著模型重數上升而顯著增加。

旋翼葉片的周期性旋轉引起回波的周期性調制，當葉片旋轉到形成鏡面反射的角度時，回波在時域形成峰值閃爍，非鏡面反射時則呈現為波谷，因此時域回波呈現為間隔一定時間出現的譜峰。在時頻域內，該閃爍表現為占據一定帶寬的頻率帶，并且頻率帶的相對位置關系與旋翼結構的轉動參數密切相關。目前常用的參數提取算法主要是根據譜峰間隔與葉片數量和轉速的關系，以及最大微多普勒頻率與葉片尺寸與轉速的關系實現。然而，本文考慮從閃爍位置變化的角度進行旋翼參數估計，然后在時頻圖譜中直接得到旋翼葉片數量。針對外輻射源旋翼參數估計問題，本文首先分析時頻域閃爍分布的數學形成機理，然后根據回波相位項導致時頻域閃爍呈現起伏的特征，提出利用相位補償的方法完成同一葉片閃爍的零頻聚焦。最終，在時頻圖譜中直接獲得旋翼目標包含的葉片數，然后根據閃爍中心頻率距離基準頻率最近的原則估計其他參數。另外，本方法利用最佳相位補償之后的時頻圖譜還可以實現旋翼目標在回波域的葉片分離。

2 GNSS外輻射源雷達旋翼目標微動模型

2.1 外輻射源雷達旋翼目標探測模型

圖1(a)所示為外輻射源雷達的目標探測模型，外輻射源雷達通常設置有兩個接收通道，觀測通道接收目標反射回波，參考通道則主要用于接收直達波。針對GNSS信號為連續波(周期Tp1 ms，碼片速率為1.023 Mbits/s)的特點，此處考慮采集時長為Tp的信號片段構成一個脈沖，然后對采集的多個脈沖進行相干處理。圖1(a)中R0,Rt和Rr分別表示衛星到接收機的直達距離、雙基地雷達發射距離和接收距離，雙基地角為β，波長為λ。設旋翼目標的葉片旋轉速率為fr，相應的旋轉角速度為w，且與角平分線夾角為φ，如圖1(b)所示。則接收回波與直達波信號匹配濾波之后可表示為

圖1 外輻射源雷達旋翼目標微動模型Fig.1 Micro-motion model of rotor targets in PBR

其中，σ0為回波幅度，Rbis(m)為雙站距離。pr(n,m)表示匹配濾波之后包絡函數。

2.2 旋翼目標微動模型

本文設定外輻射源地雜波和機身回波已利用經典算法進行抑制，例如擴展相消算法[1]、自適應復變分模態分解[18]、經驗模態分解[19]等，回波信號只剩下目標旋翼回波以及噪聲。以旋翼直升機為例進行討論，并考慮目標處于懸停狀態，式(1)中雙站距離Rbis(m)可建模為

其中，Rbis0為旋轉中心的雙基地距離，l0表示散射點距旋轉中心的距離，葉片總長度為L，該葉片初始角為θ0。cos(β/2)cosφ可稱為雙基地因子，主要由目標與雷達系統的相對位置決定，此處設定已完成目標定位，并且已完成目標平動補償。

根據文獻[12]可推導外輻射源條件下單葉片的回波模型為

其中，A(n,m)表示接收信號包絡幅度及與微動無關的平動相位。

設旋翼目標共有Q個葉片，則該旋翼目標的回波信號可建模為

那么，第q個葉片旋轉運動引起的微多普勒頻率為

式(5)表明目標回波信號微多普勒頻率在慢時間域受正弦函數sin(wm+θq)的調制，即微多普勒頻率是隨慢時間呈非線性變化的。另外，式(5)表明微多普勒頻率的最大值由目標旋翼的長度、轉速、波長以及目標所處的雙基地位置共同決定。當其中相位項φwm+θq滿足式(6)時，多普勒頻率呈現最大值，且最大值如式(7)所示，為該旋翼目標的最大微動多普勒頻率。另外，結合式(4)和辛克函數的特點，此時辛克函數出現極值，即在慢時間域出現閃爍。

3 基于相位補償的旋翼目標微動特征提取

3.1 旋翼目標回波時頻分析

仔細考察式(4)，忽略回波振幅的影響，單個葉片的回波主要由辛克函數與相位的乘積構成。不失一般性，第1個葉片的回波可簡寫為

由式(5)和式(8)可見，單個葉片的回波呈現為多個辛克函數求和的形式，且回波多普勒頻率是時變的，所以此處考慮采用短時傅里葉變換(Short Time Fourier Transform,STFT)在慢時間維對回波的多普勒頻率進行分析。由傅里葉變換的性質可知：

其中，FT{ }表示傅里葉變換。鑒于sinc(t)與rect(f)為一對傅里葉變換對，且指數函數的傅里葉變換易于求解，具體如式(10)和式(11)所示：

可見，回波包絡的傅里葉變換為占據一定帶寬的矩形條帶函數，在時頻圖譜中為閃爍條帶，如圖2(a)所示；而相位的傅里葉變換為一段正弦函數，如圖2(b)所示。

由式(9)可知，單葉片回波的傅里葉變換為包絡的傅里葉變換與相位傅里葉變換的卷積，進而可得完整回波信號的時頻圖譜如圖2(c)所示。這就是旋翼葉片回波在時頻域呈現隨慢時間起伏的閃爍條帶的數學形成機理。因此，如果將該葉片回波中的相位補償完全，對應葉片的時頻圖譜將由圖2(c)變成圖2(a)，即完成該葉片的相位補償，這樣便可在多個葉片產生的多次閃爍當中將完成相位補償的葉片識別并提取出來，且聚焦到零頻的閃爍屬于相同的葉片，根據兩個零頻閃爍間隔的閃爍數量估計整個旋翼的葉片數。

3.2 相位補償算法

根據式(8)，建立回波的補償相位：

其中，α,μ,χ分別表示幅度、角速度與初相的系數向量，這樣就將旋翼參數估計問題轉換成最優化求解問題了。根據3.1節所述，單葉片回波經過相位補償之后的時頻分布為一系列中心頻率為零的零頻閃爍，下面詳細分析不同旋翼結構下的基準頻率。

3.2.1 非對稱結構相位補償

當目標旋翼包含有奇數葉片時，由于葉片的不對稱性，時頻域正負頻率閃爍交替出現，其時頻分布類似圖2(c)。鑒于回波相位的參數Λ與的關系Λ即參數α與μ相關聯。設fCD表示單葉片閃爍條帶的帶寬，則fCD因此參數搜索維度就由三維變換到二維 (μ,χ)。進一步推導，可將式(12)重寫為

圖2 單個葉片的時頻圖譜Fig.2 The TFD of signal blade

仔細考察相位項發現，其正相關于兩個正弦函數的差值，并且差值的變化趨勢與參數 (μ,χ)相關。具體分以下3種情況討論：

(1) 當wμ,θ1χ時，該葉片回波相位被補償完全，此時：

(2) 當wμ,θ1?χ時，剩余相位為兩個同頻但初相不等的正弦信號的差值，變化趨勢與正弦信號相同：

(3) 當w ?μ,θ1χ時，剩余相位為兩個不同頻的正弦信號之間的差值，在時間域呈現為幅度逐漸變化的振蕩：

為充分說明上述不同情況下的回波多普勒頻率的變化情況，分別設置不同參數，計算剩余相位exp{j2π(cos(wm+θ1)?cos(μm+χ))}的多普勒頻率，其波形如圖3所示。

圖3(a)中黑色虛線是wμ4.8×2π,θ1?χ=π/9時的結果，可見當正弦函數角頻率w相同時，兩個正弦函數的差值依然是正弦函數，因此剩余相位的多普勒頻率依然被正弦函數調制，時頻圖譜中的閃爍也無法聚焦到零頻。圖3(b)中紅線表示相位補償完全情況下的多普勒頻率，可見多普勒頻率恒等于零。黑色虛線表示w4.8×2π,μ4×2π,θ1χ時的計算結果，可見不同頻正弦函數的相位補償將導致多普勒頻率隨時間變化，且變化趨勢為幅度變化的振蕩。因此，不同頻的補償相位無法將回波相位補償完全，時頻圖譜上相同葉片的全部閃爍無法聚焦到零多普勒頻率。

圖3 不同正弦函數差值導數波形圖Fig.3 The result of different sinusoidal function

綜上，第1種情況多普勒頻率始終為零頻，第2種情況多普勒頻率呈現為規律的正弦變化曲線，第3種情況下同樣未形成聚焦，表現為幅度不斷變化的不規則振蕩。因此，僅在相位補償完全的條件下，非對稱結構旋翼目標閃爍的基準頻率為零。

3.2.2 對稱結構相位補償

當目標旋翼包含有偶數葉片時，由于葉片的對稱性，時頻域內同一時間會出現正頻率閃爍和負頻率閃爍。其中，正頻率閃爍表示閃爍條帶中心頻率為正值，負頻率閃爍的中心頻率則為負值。

進一步，兩個對稱葉片的回波可簡寫為

那么，相位補償之后的回波可表示為

仔細考察式(19)可以發現，由于兩個葉片的初始相位不同，無法同時完成兩個葉片的相位補償，所以對稱結構旋翼回波閃爍無法聚焦到零頻。當w ?μ或θ1?χ時，相位補償情況如上一節所述，因此本部分主要討論最佳補償情況。當wμ,θ1χ時，式(19)可推導為

由式(20)可知，在產生閃爍的時間，兩個對稱的葉片分別產生零頻閃爍和非零頻閃爍。具體后者的多普勒頻率如式(21)所示。

根據式(21)可知，后者的多普勒頻率呈現正弦變化的趨勢，因此后者產生閃爍的中心頻率也成正弦函數起伏。考慮到該葉片對的閃爍恰好出現在同一時刻，在時頻圖譜上表現為同一時刻連接在一起的雙閃爍。其中一個葉片的相位被補償完全，該葉片的回波只剩慢時間域包絡，時頻圖譜中的閃爍被聚焦到零頻附近。另一個葉片的相位變為原來的負二次方，其聚焦的多普勒頻率變為。顯然，該葉片對形成的閃爍連接在一起構成雙閃爍，且該雙閃爍的中心頻率位于處。

總體而言，該時刻雙閃爍的頻率帶寬為單閃爍的兩倍，并且中心頻率由零變成振幅為Λw的正弦函數。當該對稱葉片旋轉半圈之后，即相位增加π，閃爍的中心頻率則由?Λw變化到Λw，然后再次變化到?Λw。雖然此時閃爍的中心頻率隨時間變化不定，但同一對稱葉片產生雙閃爍的中心頻率僅僅出現Λw和?Λw兩種結果。因此仍然可以將時頻域閃爍條帶的中心頻率聚焦到確定值，此處選擇基準頻率為負值?Λw。

3.3 參數提取最優化準則

針對上述分析，為克服剩余相位多普勒頻率的不規則振蕩并減弱單個閃爍聚焦效果的不穩定性引起的估計誤差，此處考慮利用回波時間內單葉片或者對稱葉片對的所有閃爍進行最優化參數搜索。另外，此處假設同一旋翼目標的所有葉片轉速相同，差別僅在于各葉片的初相。因此，獲得其中一個葉片參數即可獲得整個旋翼的相關參數，所以本算法考慮利用單個葉片或者葉片對的相位補償實現整個旋翼目標的特征提取。由于相位補償之初無法得知哪些閃爍屬于同一葉片，并且可能出現多個補償相位使某個閃爍聚焦到相同的多普勒頻率，所以此處考慮采用迭代的方式逐漸搜索到屬于某一葉片的全部閃爍，并將其聚焦到零頻，并以此為原則估計目標的各項參數。

當旋翼結構非對稱時，基準頻率為零，此時計算閃爍中心頻率的絕對值即表示閃爍中心距離基準頻率的頻率差。起始階段可以選擇時頻圖譜中的任意一次閃爍作為基準閃爍([D1])，利用基準閃爍中心頻率最小進行第1次迭代，搜索獲得最佳補償之后的時頻圖譜。然后再利用閃爍中心頻率最小原則對基準閃爍進行擴展，如3.2.1節所分析。根據其余各閃爍中心頻率最小的原則尋找第2個閃爍([D2])，并加入基準閃爍([D1,D2])，然后根據該兩個閃爍進行第2次迭代，并完成基準閃爍的更新；如此逐漸增加基準閃爍([D1,D2,...,DI])的數量，直至剩余其他閃爍的中心頻率均大于閾值。此時觀測時間內該葉片所有閃爍完成聚焦，而其他葉片的閃爍未完成聚焦。最佳迭代之后的相位補償結果在時頻圖譜上呈現為一系列規律的零頻閃爍，然后根據兩次零頻閃爍之間間隔的閃爍數確定旋翼目標包含的葉片數量，并估計轉速和初相。

根據式(7)可知，在獲得目標的轉速之后，可以根據閃爍條帶的多普勒帶寬估計旋翼葉片的長度。

同樣，當旋翼結構對稱時，對稱葉片產生的雙閃爍基準頻率為?Λw，根據上述逐步聚焦的原則，將所有應當聚焦到?Λw的雙閃爍聚焦完成。此外，如3.2.2節所述，兩個聚焦到?Λw處的雙閃爍之間的位于Λw處的雙閃爍也屬于該葉片對。然后從時頻圖譜中直接獲得旋翼葉片數，并根據整體雙閃爍中心頻率距離?Λw最近的原則估計該旋翼目標的其他相關參數。

綜上，通過相位補償算法，可以實現旋翼目標的特征參數提取，具體算法執行框圖如圖4所示，其主要步驟如下：

圖4 算法執行框圖Fig.4 The flowchart of the proposed method

(1) 采集目標原始數據，并進行距離向壓縮獲取目標所在距離單元的慢時間回波；

(2) 對原始數據進行時頻分析，確定葉片結構對稱與否；

(3) 根據葉片是否對稱分別確定基準頻率；

(4) 任意選取時頻圖譜中的一次閃爍作為基準閃爍；

(5) 根據設置的參數矩陣(μ,χ)構造補償相位依次對慢時間回波進行相位補償；

(6) 獲取每次補償回波的時頻圖譜，根據基準閃爍的中心頻率距離基準頻率最近的原則確定當前最優時頻圖譜；

(7) 遍歷步驟(6)獲得的最佳時頻圖譜中剩余所有閃爍，選取中心頻率距離基準頻率最近的閃爍加入基準閃爍，并回到步驟(5)；

(8) 直至時頻圖譜中所有閃爍的中心頻率距離基準頻率大于閾值，迭代循環結束，基準閃爍不再更新；

4 仿真實驗

本節結合阿帕奇AH-64武裝直升機旋翼[12]模型，利用GPS L1波段仿真信號進行實驗，具體仿真參數如表1所示。為充分說明所提參數估計方法的有效性和對于不同參數空間設置的適應性，本文所提方法的仿真實驗效果還與其他經典方法進行了比較。

4.1 基于相位補償的參數提取方法實驗結果

4.1.1 非對稱旋翼目標實驗

本部分仿真實驗設置目標旋翼包含3個葉片，其他參數如表1所示。圖5為旋翼目標的回波及時頻分布，其中圖5(a)為目標所在距離單元的慢時間域回波。可見，因為目標葉片的不斷旋轉，慢時間域回波主要由一系列相同間隔的譜峰構成，稱之為閃爍。圖5(b)為慢時間域回波序列的短時傅里葉變換，顯然正頻率閃爍和負頻率閃爍交替出現在不同時間。另外，可以根據圖5(b)估計得到單次閃爍的帶寬為284.099 Hz。由前述分析可知，閃爍在時頻域的不對稱性主要由旋翼結構的不對稱性引起，因此可以確定目標旋翼包含奇數葉片。

圖5 非對稱旋翼目標回波分析Fig.5 The echoes of asymmetry rotor targets

表1 仿真實驗參數Tab.1 The parameters in simulation experiments

為確定目標旋翼包含的葉片數，采用相位補償的方式對目標慢時間域回波進行處理，然后再進行時頻分析，圖6為相應的時頻分析結果。本實驗選擇慢時間域上第2個閃爍(0.067 s)作為原始基準閃爍，然后根據基準閃爍的中心頻率距離基準頻率最近的原則選取最佳補償的時頻圖譜，其中圖6(a)為利用原始基準閃爍進行相位補償獲得的最佳聚焦效果。由圖6可見，原始基準閃爍的中心頻率接近零，而圖6中的其他閃爍則在一定程度上偏離了零頻，即僅僅原始基準閃爍被聚焦到零頻。如前3.2.1節所述，此現象是由補償相位的角頻率不同于葉片旋轉的角頻率導致，相位補償之后的多普勒頻率在慢時間域呈現為不同幅度的振蕩，與圖3呈現的結果相吻合。為克服單個閃爍對聚焦效果穩定性的影響，本方法考慮利用多個閃爍進行聯合聚焦。圖6(b)所示為第2次迭代獲得的最佳聚焦效果，此時基準閃爍包含有兩個不同位置接近于零頻的閃爍(0.067 s,0.171 s)。觀察圖6(b)發現，前3個閃爍都被聚焦于接近零頻的位置，而其他閃爍則明顯偏離零頻。圖6(c)和圖6(d)分別表示第3、第4次迭代獲得時頻圖譜，其基準閃爍分別包括3個閃爍(0.067 s,0.171 s,0.275 s)和4個閃爍(0.067 s,0.171 s,0.275 s,0.379 s)。迭代過程中基準閃爍的更新情況如表2所示。觀察發現兩圖中各閃爍的分布結構大致相同，白色箭頭所指的5個不同位置的閃爍均被聚焦于零頻，因此該5個閃爍(0.067 s,0.171 s,0.275 s,0.379 s,0.483 s)可以確定為同一葉片在不同時間所產生。結合圖6(c)和圖6(d)可以發現，相同葉片兩次不同時間的閃爍之間還間隔兩次其他葉片的閃爍，這樣就可以判斷該旋翼目標包含有3個葉片。

表2 非對稱旋翼結構基準閃爍的更新Tab.2 The update of standard flashes in asymmetry rotor target

為更直觀呈現相位補償過程中前述5個閃爍中心頻率的變化，圖7給出了4次迭代過程中該5個閃爍中心頻率的變化。其中圓圈狀藍色虛線表示第1次迭代之后該5個閃爍中心頻率的變化，可見不同位置閃爍的中心頻率變化較大，僅基準閃爍距離基準頻率比較近；菱形紫色線表示第2次迭代之后該5個閃爍中心頻率的變化，相比前一次，中心頻率起伏變小，同樣僅有基準閃爍距基準頻率較近；依次矩形黑色虛線和三角形紅色線分別表示第3、第4次迭代之后該5個閃爍中心頻率的變化，兩次迭代之間曲線變化不太明顯，且基本等于零。可見，隨著基準閃爍的逐漸更新，基準閃爍包含的閃爍數量也逐漸增加。圖7表明隨著基準閃爍數量的增加，相同葉片的閃爍越接近于基準閃爍。圖8為最后一次迭代之后的參數搜索結果，其中Z軸表示基準閃爍的中心頻率距離基準頻率的平均值，X軸表示的是初始相位，Y軸表示葉片轉速。圖6(d)所示的最佳相位補償效果對應的參數如圖8中指數所示，可得到目標旋翼的轉速為4.80 r/s。另外，根據式(22)得目標旋翼葉片長度為7.34 m。

圖8 最佳參數搜索結果Fig.8 The optimal parameters search result

4.1.2 對稱旋翼目標實驗

前述驗證的是非對稱結構旋翼目標的特征提取效果，由于對稱結構旋翼產生的閃爍與非對稱結構存在一定的差異，本部分用于驗證所提方法對于對稱結構旋翼目標的特征提取效果，實驗目標為相同葉片長度、轉速的旋翼，但包含4個葉片。

圖9(a)所示為4個葉片旋翼目標原始回波的時頻圖譜，相比于圖5(b)而言，此處時頻圖譜中的閃爍是對稱結構，同一時刻的閃爍既包括正頻率，又包括負頻率，稱為雙閃爍，原因如3.2.2節所述。根據圖9(a)可估計雙閃爍所占帶寬為571.73 Hz，此時fCD285.87 Hz。此處同樣選擇第2個雙閃爍(0.067 s)作為原始基準閃爍，圖9(b)為利用原始基準閃爍獲得的最佳補償之后的時頻圖，可估計其基準閃爍的中心頻率為—141.28 Hz，約為雙閃爍帶寬的1/4，仔細觀察時頻圖可見，該圖中的雙閃爍并不存在明顯的規律性。圖9(c)和圖9(d)分別表示第2、第3次迭代之后的時頻圖，過程中所用的基準閃爍如表3所示。兩次聚焦效果基本相同，表明此處只要結合兩個雙閃爍進行參數搜索，就可以實現比較精確的相位補償。仔細觀察時頻圖9(c)可以發現，經過相位補償之后，首先可以確定第2次(0.067 s)、第6次(0.275 s)、第10次(0.483 s)雙閃爍屬于同一葉片對，如圖9中白色箭頭所示。其次考慮到同一葉片對所形成雙閃爍的中心頻率僅呈現兩種對稱的結果，第4次(0.171 s)、第8次(0.379 s)雙閃爍的中心頻率為第2次雙閃爍的中心頻率的負值，可見第4次、第8次雙閃爍也與第2次雙閃爍同屬同一葉片對。因此，分析可知圖9(d)中的第2次、第4次、第6次、第8次、第10次雙閃爍屬于同一葉片對，具體如圖9中白色箭頭所指示。此外，同一葉片對產生的雙閃爍之間僅間隔一個雙閃爍，且考慮到每個雙閃爍對應兩個葉片，那么可以獲得該旋翼結構包含4個葉片。

圖9 對稱旋翼結構相位補償之后的時頻圖譜Fig.9 The TFD of symmetry rotor blade after phase compensation

表3 對稱旋翼結構基準閃爍的更新Tab.3 The update of standard flashes in symmetry rotor target

同樣，圖10給出了3次迭代過程中基準閃爍中心頻率的變化，其中紅色實線表示相位補償的基準頻率。圖中圓圈狀藍色虛線表示第1次迭代之后中心頻率分布情況，由于第1次選擇的是0.067 s處的基準閃爍，所以僅僅使得該雙閃爍的中心頻率接近于基準頻率，而其他雙閃爍的中心頻率則距基準頻率較遠。另外，圖10中菱形紫色線表示第2次迭代之后中心頻率的分布，其3個雙閃爍的中心頻率都相同程度地非常接近于基準頻率。黑色矩形虛線表示第3次迭代之后中心頻率的分布情況，其結果與第2次迭代幾乎重合，可見在第2次迭代之后即可獲得最佳相位補償效果，即利用同一葉片兩個不同位置的雙閃爍作為迭代最優化準則基本可獲得最佳的參數估計結果。

圖10 迭代過程中基準閃爍中心頻率變化Fig.10 The variety of center frequency of standard flashes during iterations

圖11為經過最后一次迭代之后的參數搜索結果，由圖11可見，基準閃爍中心頻率距離基準頻率最小值為2.00 Hz，幾乎等于多普勒分辨率(2 Hz)，且達到最小值時葉片轉速為4.80 r/s，相位為5.93 rad。因為不同葉片之間僅初始相位不同，所以此處關注的參數為葉片轉速，與表1中參數設置相比，可見參數估計精度較高。然后綜合目標所在的雙基地參數，根據式(22)可計算獲得葉片的長度為7.38 m，可見參數估計精度也比較高。

圖11 最佳迭代之后的參數搜索結果Fig.11 The optimal parameters search result

4.1.3 參數估計精度分析

上述部分已經驗證了本算法對于對稱旋翼結構和非對稱旋翼結構特征提取的有效性，下面進一步分析所提算法的參數估計精度。旋翼葉片數的估計誤差直接決定算法的有效性，而且本算法利用時頻圖譜中基準閃爍之間的位置關系實現參數估計。另外，3.2.1節已經分析過補償相位的初始相位差僅會帶來閃爍的整體平移，而不會影響閃爍之間的相對位置關系。葉片長度的估計一方面受閃爍條帶帶寬估計的影響，另一方面也受旋轉速度的影響，所以此處主要討論旋翼轉速的估計精度。圖12給出了不同轉速相位補償之后獲得的最佳時頻圖譜中基準閃爍的中心頻率距離基準頻率的差值，然后重復50次蒙特卡洛實驗取平均值。圖12中紅色矩形實線轉速為4.8 r/s的結果，其5個位置閃爍的頻率差明顯低于其他轉速條件下的頻率差，所以當補償相位的轉速與實際轉速相等時可獲得最佳聚焦效果。并且距離紅色矩形實線最近的分別是轉速為4.7 r/s和4.9 r/s的聚焦效果，其次是轉速是4.6 r/s和5.0 r/s的聚焦效果。可見，當補償相位的轉速越接近目標實際轉速時，其獲得的補償效果越好，因此利用本文所提出的相位補償算法可以實現比較精確的旋翼目標特征提取。

圖12 不同角頻率補償相位對閃爍位置的影響Fig.12 The effect of different angular velocity on locations of standard flashes

4.1.4 葉片分離效果

本文所提算法通過相位補償的方式將某個葉片的相位補償完全，然后該葉片的閃爍就聚焦在零多普勒頻率處。然后根據時頻圖譜直接獲得旋翼包含的葉片數，并估計其他相關參數。此外，該算法還能實現旋翼葉片或葉片對的分離。以三葉片旋翼目標為例，圖13為利用相位補償之后的時頻圖譜進行葉片分析的效果。以相位補償完全的葉片為基準，將其視為第1個葉片，如圖13(a)中紅色箭頭所示；顯然，根據各閃爍距離零頻閃爍的時間維距離，依次可以確定第2個葉片(白色箭頭所示)和第3個葉片(黃色虛線箭頭所示)。對應到慢時間域回波，其時域閃爍可分離為相互不同的3個葉片形成，相同葉片形成的閃爍之前間隔兩個其他葉片的閃爍，如圖13(b)所示。

圖13 葉片分離的效果Fig.13 The blade separation effect

為進一步說明所提方法對旋翼葉片分離的效果，圖14和圖15分別為對圖5(a)和圖6(d)進行葉片分離的結果，具體是利用峰值提取技術將各葉片形成的閃爍分別提取出來。前述部分已經估計得到旋翼包含3個葉片，所以分片分離后主要包含3個葉片，各葉片的起始相位自左至右分別為：1.75π,0.42π,1.09π，每個葉片的慢時域回波如圖14所示。由圖中可得各葉片譜峰之間的時間間隔約為0.105 s，為圖5(a)中所示譜峰間隔(0.035 s)的3倍，與旋翼包含三葉片的結果相符。圖15為各葉片回波的時頻圖譜，可見每個葉片時頻圖譜都呈現為正頻率閃爍和負頻率閃爍交替出現的情況。雖然形成閃爍均表示葉片與雷達系統之間形成鏡面反射，但是該幾何關系的形成可以分為兩種情況：葉片靠近雷達運動和背離雷達運動。具體葉片運動方法與正負閃爍的對應關系與發射信號的相關，此處不作詳細討論。所以，本文所提方法不僅能實現旋翼參數的精確估計，而且可以實現不同葉片回波的分離，且能對應到不同的初始相位。

圖14 葉片回波分離Fig.14 The echoes of different blade

4.2 其他經典算法效果對比

此處設置與4.1.1節相同的實驗環境，設置旋翼包含3個葉片，利用經典的特征提取方法對信號進行處理。圖16表示對慢時間域回波直接進行傅里葉變換處理的結果，可見回波頻譜大致分布在—273.5～273.5 Hz帶寬的范圍內，與4.1.1節估計結果相近，但是該結果僅給出了多普勒頻率范圍，而沒有呈現多普勒頻率的變化趨勢，因此難以提取出回波中的微多普勒頻率信息。圖17為Hough變換對原始時頻圖譜進行正弦曲線檢測的結果，其中X軸表示轉速，Y軸表示初始角度。可見轉速估計結果依次為(4.9 r/s,4.7 r/s,4.7 r/s)，均值為4.77 r/s，進一步計算得到葉片長度為7.36 m。可見，受時頻分析工具STFT分辨率的限制，Hough變換參數空間結果中的3個峰值都不夠銳利，容易影響峰值檢測的精度，進而導致參數估計結果精度比較低。

圖15 時頻圖中的葉片分離Fig.15 The TFD of echoes of different blade

圖16 傅里葉變換的結果Fig.16 The result of Fourier transform

圖17 Hough變換的結果Fig.17 The result of Hough method

圖18為不同參數設置條件下的OMP算法實驗結果，具體不同實驗的參數設置如表4所示，圖中X軸表示葉片長度，Y軸表示目標葉片的轉速，Z軸則表示幅度。其中，圖15(a)對應實驗a的結果，可見圖中包含3個峰值，可取該3個峰值指數的平均值作為實驗a的參數估計結果。由于3個峰值分別位于不同的旋轉速度上，即旋轉速度估計結果不一致，所以對應的葉片數也就不相等，可認為該條件下的OMP算法不能獲得目標的參數估計。雖然此時仍然可以通過3個峰值的均值計算出參數估計結果，但是均值計算結果沒有意義，此處將無效的葉片長度與轉速的估計結果寫在表5內。圖18(b)表示實驗b的結果，由圖中指數可知，3個峰值對應相同的轉速4.80 r/s，對應葉片數為3。此時可以認為本條件下的OMP參數估計方法可以有效估計目標參數，進而可得目標旋翼的長度估計結果為7.45 m，具體實驗b的有效參數估計結果如表5所示。圖18(c)為實驗c的結果，此處設置初相角在0～120°內分10等份。與圖18(a)類似，此處3個峰值也對應不同的旋轉速度，同樣可以認為實驗c參數設置條件下OMP算法無法獲得有效的目標參數估計結果，實驗c的相應無效結果也寫在表5中。根據表5可知，當OMP算法無效時，其參數估計結果遠遠偏離真實值。因此，為保證OMP算法的有效性，該算法需要設定合適參數搜索區間和網格劃分疏密。

圖18 OMP參數估計方法的實驗結果Fig.18 The experiment result of OMP algorithm

表4 OMP算法實驗條件Tab.4 The experiment condition of OMP algorithm

表5 參數估計結果對比Tab.5 The parameter estimation result of different algorithms

對比實驗a與實驗b的參數設置，發現基于OMP的參數搜索算法容易收斂到局域的極值，而不能穩定的收斂到全域的最值，導致參數估計結果可靠性較差。同樣，比較實驗b與實驗c的參數設置，發現參數空間的網格化程度也會嚴重影響到OMP方法參數估計結果的有效性。

表5給出了不同實驗條件下的參數估計結果。本文所提方法無論旋翼結構是否對稱都能夠實現有效的參數估計，且相比于Hough變換和OMP算法的參數估計精度更高。另外，基于OMP算法的參數估計僅能在實驗b的參數設置條件下獲得有效的估計結果，其他兩個實驗條件下無法獲得有效的估計結果。綜上所述，本文所提出的相位補償算法相比Hough變換和OMP算法而言，可以實現性能更加穩定且結果更加準確的參數估計。另外，相比于其他經典算法而言，該算法還可以實現旋翼目標在回波域的葉片分離，為后續進一步的旋翼目標識別提供更多思路。

5 結論

本文分析了外輻射源旋翼目標微動信號模型，并結合時頻域閃爍特征和微動信號特性，提出了一種基于相位補償的旋翼目標參數估計方法。首先利用時頻分析方法獲得回波的時頻圖譜，然后根據閃爍的形成機理對回波信號進行相位補償，并根據相位補償結果估計得到目標的葉片數量、長度以及轉速。在算法執行過程中，為獲得更高的算法穩定性，通過迭代的方式逐漸增加閃爍的數量進行聯合搜索，避免單個閃爍的不穩定性對算法性能的影響。另外，針對對稱結構和非對稱結構旋翼的不同特點，分別采用不同的基準頻率，使得本方法對不同旋翼結構的適應性更強。最后利用仿真實驗驗證了本文所提方法具有更高的可靠性和更準確的參數估計結果，并且該方法可以有效實現旋翼目標在回波域的葉片分離，為后續進一步的旋翼目標識別提供更多的解決思路。