基于線性調(diào)頻信號二次相位變換的主動聲吶目標(biāo)速度估計

顧怡鳴，宮在曉，李整林

(1. 中國科學(xué)院聲學(xué)研究所聲場聲信息國家重點實驗室，北京 100190；2. 中國科學(xué)院大學(xué)，北京 100190)

0 引 言

線性調(diào)頻(Linear Frequency Modulation, LFM)是主動聲吶常用的探測信號，通過分析 LFM 信號目標(biāo)回波可估計目標(biāo)與探測平臺的相對速度。傳統(tǒng)的LFM測速方法主要有多普勒頻移補償法和組合脈沖法。多普勒頻移補償法利用長脈寬 LFM 信號的多普勒敏感性，對探測信號進(jìn)行多普勒頻移搜索來估計相對速度[1]。組合脈沖法利用LFM信號模糊函數(shù)的時延-多普勒耦合性，通過正負(fù)調(diào)頻信號匹配濾波峰值的時延差估計速度，也被擴展應(yīng)用于雙曲調(diào)頻信號[2]。但主動聲吶的發(fā)射占空比通常較小，組合脈沖增加了信號的發(fā)射周期。

分?jǐn)?shù)階傅里葉變換(Fractional Fourier Transform, FrFT)是一種時頻旋轉(zhuǎn)變換，近年來已經(jīng)被廣泛用于 LFM 信號的檢測和參數(shù)估計[3-13]。文獻(xiàn)[14]首先將 FrFT應(yīng)用于主動聲吶回波檢測，并利用海上實驗數(shù)據(jù)證明了 FrFT在混響背景下可以獲得比匹配濾波器更高的信噪比和速度估計精度。但FrFT較高的計算復(fù)雜度限制了其實際應(yīng)用，針對該問題，學(xué)者們提出了一些快速 FrFT算法[15-21]。采樣型算法對時間和頻率進(jìn)行尺度歸一化，可以借助快速傅里葉變換(FFT)算法極大地降低計算量[15-16]；特征向量分解型算法[12-13]通過設(shè)計離散傅里葉變換矩陣求解對應(yīng)的特征向量；文獻(xiàn)[21]中利用莫比烏斯反演公式證明了分?jǐn)?shù)階傅里葉級數(shù)的存在，并以此為基礎(chǔ)提出了低復(fù)雜度的離散FrFT算法。FrFT方法檢測 LFM 信號主要利用了其解調(diào)頻的性質(zhì)，舍棄一些冗余的性質(zhì)，對變換形式進(jìn)行簡化以降低計算量，例如二次相位變換[22-24](Quadratic Phase Transform, QPT)或簡明分?jǐn)?shù)階傅里葉變換[7-8]，以及具有高次變換核的多項式時頻變換[25-26]。

現(xiàn)有的 FrFT估計速度方法需要對變換參數(shù)進(jìn)行搜索，在自主水下機器人(Autonomous Underwater Vehicle, AUV)等計算力有限的平臺難以實時處理。本文針對該問題提出了一種基于二次相位變換的線性調(diào)頻信號速度估計方法。QPT參數(shù)失配時能量無法聚焦，幅度譜存在一定展寬。本文推導(dǎo)了 QPT參數(shù)與信號參數(shù)失配時變換域幅度譜的解析解，分析了二次相位變換幅度譜展寬與多普勒頻移之間的關(guān)系，提出了利用單次QPT后幅度譜寬度的速度估計方法?；趩蝹€線性調(diào)頻信號進(jìn)行QPT后幅度譜的展寬信息實現(xiàn)了目標(biāo)相對速度估計，避免了對變換參數(shù)的掃描，減少了一個維度的計算量。針對單次QPT結(jié)果無法確定速度符號的問題，提出了對稱參數(shù)QPT符號判斷方法。通過數(shù)值仿真驗證了方法的有效性，并與 FrFT參數(shù)掃描方法進(jìn)行了對比。

1 FrFT測速方法

FrFT的變換參數(shù)cotα與信號的調(diào)頻斜率k匹配時，變換核中的 ( t2/2)cotα項對信號實現(xiàn)解調(diào)頻，變換結(jié)果的幅度取得最大值。文獻(xiàn)[14]依據(jù)該性質(zhì)給出了通過掃描 FrFT變換參數(shù)搜索最大峰值的測速方法。以長度為20 s、頻帶400～500 Hz的LFM信號為例，根據(jù)式(5)得到FrFT對LFM信號u-多普勒模糊(見圖1)。由圖1可以看出，當(dāng)變換參數(shù)與信號速度匹配時變換結(jié)果的幅度取得最大值。FrFT測速方法的原理是在變換域-多普勒平面內(nèi)搜索最大值，找到最大值對應(yīng)的變換參數(shù)，進(jìn)而通過變換參數(shù)與多普勒頻移之間的關(guān)系得到目標(biāo)相對速度的估計值。

圖1 FrFT對LFM信號的u-多普勒模糊圖Fig.1 Schematic diagram of the u-Doppler ambiguity of FrFT for LFM signal

2 QPT速度估計方法

2.1 參數(shù)失配時LFM信號的QPT

二次相位變換可以看作是 FrFT的簡化，信號x(t)的QPT定義為[25]

其中，u和κ是變換參數(shù)。當(dāng)變換參數(shù)和線性調(diào)頻信號的調(diào)頻斜率k匹配，即k =κ 時，信號能量會在變換域聚集為峰值，此類情況已經(jīng)在現(xiàn)有文獻(xiàn)中廣泛討論；當(dāng)變換參數(shù)與信號參數(shù)失配時k≠κ ，信號的能量在變換域表現(xiàn)為一定展寬，已有文獻(xiàn)尚未發(fā)現(xiàn)有關(guān)變換參數(shù)與信號參數(shù)失配時特征的研究。以下給出變換參數(shù)與 LFM 信號調(diào)頻斜率不匹配情況下變換結(jié)果的解析解。

多普勒頻移后的回波信號二次相位變換可表示為

變換后信號的幅度譜由兩種形式的 Fresnel積分(如圖2所示)決定，式(10)對應(yīng)的幅度譜大部分區(qū)域為零，如圖3所示，能量集中在一定的寬度，定義幅度譜寬度Bα為幅度譜峰值-3 dB對應(yīng)的u域?qū)挾?，能量集中的范圍的寬度?/p>

圖2 宗量為x的Fresnel積分Fig.2 Fresnel integral with argument x

圖3 多普勒頻移后的回波信號QPT的幅度譜示意圖Fig.3 QPT amplitude spectrum of the echo with Doppler frequency shift

2.2 單次QPT速度估計方法

以上推導(dǎo)表明，參數(shù)失配時LFM信號QPT的幅度譜寬度與信號多普勒頻移存在確定關(guān)系，利用該特征可估計目標(biāo)相對速度。若變換參數(shù)為原始信號參數(shù)時κ= k，LFM信號QPT的幅度譜的寬度為

式中：c為聲速；k和T為LFM信號的調(diào)頻斜率和脈沖寬度。只需要通過估計幅度譜寬度Bk即可得到相對速度。QPT幅度譜的寬度信息可通過閾值檢測獲取，即：

其中：d為檢測門限，根據(jù)閾值檢測幅度譜的前后沿u0和u1可估計幅度譜寬度的估計值目標(biāo)靜止時QPT幅度譜為sinc函數(shù)，幅度譜寬度為。需要注意的是，式(10)并不是完美的矩形窗函數(shù)，檢測門限較高時估計的寬度會偏小，需要通過修正系數(shù)進(jìn)行補償。

根據(jù) QPT幅度譜的寬度可以估計相對速度 v的大小，相對速度的符號則可以根據(jù)幅度譜前沿與的關(guān)系來確定：

假設(shè) LFM 信號長度為 20 s，頻帶 400～500 Hz，對相對速度為 0和 30 m·s-1的信號進(jìn)行QPT后得到的幅度譜分別如圖4(a)和圖4(b)所示，可以看出，當(dāng)相對速度為0時，信號能量在變換域聚焦，信號QPT幅度譜幾乎沒有展寬。相對速度為30 m·s-1時，信號QPT幅度譜展寬明顯，并且與式(11)理論分析的寬度一致。圖4(c)為QPT對LFM信號u-多普勒模糊，可以看出，隨著相對速度的增加，QPT幅度譜逐漸展寬。同時，當(dāng)相對速度為正時，幅度譜分布位于f0- kτ0右側(cè)，相對速度為負(fù)時位于f0- kτ0左側(cè)。由于信號的能量一定，QPT幅度譜的幅度隨展寬而降低。經(jīng)過閾值檢測后的寬度圖像如圖4(d)所示，根據(jù)寬度可以進(jìn)一步估計目標(biāo)相對速度。需要注意，單次QPT方法依賴式(14)判定速度的正負(fù)，實際應(yīng)用中不能準(zhǔn)確獲取時延τ0，因此單次 QPT提供的信息在判斷速度的正負(fù)方面存在缺陷。

圖4 單次QPT的速度估計方法Fig.4 Velocity estimation method with single QPT

2.3 對稱參數(shù)QPT速度符號判斷方法

單次 QPT譜寬度僅提供了相對速度的大小信息，不能確定相對速度的符號，這里給出對稱參數(shù)QPT速度符號判斷方法。圖 5為符號判斷方法的示意圖，給定一個初始寬度B0，若信號不存在多普勒頻移時，B0對應(yīng)的兩組變換參數(shù)分別為分別對回波信號進(jìn)行變換參數(shù)為κ+和κ-的 QPT，分別可以得到兩組變換的幅度譜寬度Bκ+和Bκ-，兩次對稱參數(shù) QPT幅度譜寬度為

圖5 對稱參數(shù)QPT速度符號判斷示意圖Fig.5 Schematic diagram of speed symbol judgment with symmetrical parameter QPT

由式(15)可以看出，當(dāng)信號相對速度為正時，伸縮因子η＞ 0，Bκ+＜ Bκ-；當(dāng)信號相對速度為負(fù)時，伸縮因子η＜0，Bκ+＜ Bκ-。通過對比兩次對稱參數(shù)QPT幅度譜寬度可判斷相對速度的符號。

3 仿真實驗與性能分析

本節(jié)通過數(shù)值仿真實驗驗證 QPT速度估計方法的有效性，并將QPT方法與FrFT參數(shù)掃描方法的性能進(jìn)行對比。仿真環(huán)境為無限空間，海水聲速為 1 500 m·s-1。仿真實驗的部署如圖 6(a)所示，藍(lán)圈靜止位于坐標(biāo)原點(0, 0)的收發(fā)合置的主動聲吶平臺，黑圈代表目標(biāo)起始位置(-10 km, 5 km)，虛線代表目標(biāo)自東向西以速度 200 m·s-1勻速運動的航線。根據(jù)目標(biāo)運動速度和相對平臺的方位可以得到目標(biāo)與平臺的相對速度如圖6(b)所示，正速度表示目標(biāo)接近平臺。探測信號為長度 20 s、頻帶 400～500 Hz的LFM信號，信號發(fā)射周期為60 s。在整個20 km的運動距離內(nèi)，總計發(fā)射16個脈沖。

圖6 仿真實驗部署Fig.6 Deployment of simulation experiment

仿真接收信號的采樣頻率為5 000 Hz，假設(shè)背景噪聲為高斯白噪聲，回波信號頻帶內(nèi)的信噪比為20 dB。仿真回波信號直接經(jīng)過QPT后的幅度譜如圖7(a)所示，回波QPT后幅度譜能量集中的位置與回波的時延有關(guān)。為了更加直觀地表現(xiàn)相對速度對 QPT幅度譜的影響，將 u坐標(biāo)變換為u′坐標(biāo)，結(jié)果如圖 7(b)所示。對比真實相對速度(圖 6(b))可以觀察到 QPT幅度譜寬度與相對速度大小的對應(yīng)關(guān)系，設(shè)置檢測門限為，進(jìn)行寬度檢測并依據(jù)式(12)的函數(shù)關(guān)系即可估計目標(biāo)的相對速度。圖8為對稱參數(shù)QPT得到的幅度譜寬度檢測結(jié)果，其中藍(lán)線和紅線分別對應(yīng)正參數(shù)和負(fù)參數(shù)變換對應(yīng)的幅度譜寬度，通過對比兩者之間的大小可以判斷速度的正負(fù)。

圖7 仿真數(shù)據(jù)QPT譜寬度估計Fig.7 Estimation of the of QPT spectrum width of simulation data

圖8 對稱參數(shù)QPT幅度譜寬度估計結(jié)果Fig.8 Estimation of the amplitude spectrum width by using QPT of symmetric parameters

在實際應(yīng)用中，F(xiàn)rFT測速方法不需要對0～π所有的變換參數(shù)進(jìn)行搜索，通過速度和變換參數(shù)的關(guān)系限定參數(shù)掃描的范圍可以降低計算量。速度掃描范圍確定后，計算量與掃描的分辨率相關(guān)，掃描分辨率也直接決定了 FrFT方法的速度估計精度。仿真中，F(xiàn)rFT方法的速度掃描范圍限定在-30～30 m·s-1范圍內(nèi)，掃描間隔分別采用3、1和0.3 m·s-1，對應(yīng)的FrFT點數(shù)Nα分別為10、60和100點。分別利用 FrFT和 QPT兩種方法對相對速度進(jìn)行估計，得到速度誤差如圖9所示。由圖9可以看出，F(xiàn)rFT方法的速度估計精度主要與參數(shù)掃描的分辨率有關(guān)，仿真數(shù)據(jù)單次QPT方法的速度估計精度基本可以保持在1 m·s-1以內(nèi)，與60點的FrFT方法精度相當(dāng)。

圖9 FrFT參數(shù)掃描方法與單次QPT方法速度估計誤差對比Fig.9 Comparison of velocity estimation error between FRFT parameter scanning method and QPT method

FrFT類算法的計算復(fù)雜度分析在許多文獻(xiàn)中已有討論[7-8,15,26]，盡管QPT簡化了變換核，但計算復(fù)雜度仍與 FrFT的量級相當(dāng)。忽略計算復(fù)雜度的低階項，兩種變換基于 FFT的快速算法計算復(fù)雜度量級均為O(NlgN)，其中N表示單次處理的時域數(shù)據(jù)點數(shù)。遠(yuǎn)程主動探測回波處理的時間窗較長，處理點數(shù)的增加使算法的時間復(fù)雜度急劇上升。以第2節(jié)仿真實驗中采用的參數(shù)(窗長度40 s，采樣點數(shù)20 000點)為例，評估不同參數(shù)掃描量的計算量，結(jié)果如圖10所示，Nα=1對應(yīng)的虛線表示單次QPT的計算量。如果只估計相對速度大小，QPT計算量僅為精度相當(dāng)?shù)腇rFT方法的1/60；增加判斷相對速度符號，則計算量為FrFT方法的1/30。

圖10 不同參數(shù)掃描量的計算量估計Fig.10 The computational costs of different scanning parameters

參數(shù)失配時能量在變換域并非聚焦于一點，因此 QPT方法的抗噪聲性能差于 FrFT參數(shù)掃描方法。假設(shè)回波信號的帶內(nèi)信噪比為 20 dB，不同相對速度回波信號的QPT幅度譜波形如圖11所示，相對速度為0時，信號能量完全聚焦，背景起伏很低，隨著相對速度提高，背景起伏逐漸升高。定義峰值信噪比為輸出信噪比，通過 20次蒙特卡洛仿真估計QPT和FrFT兩種方法的輸出信噪比隨相對速度的變化，結(jié)果如圖12所示，藍(lán)線和紅線分別表示QPT方法和FrFT方法的輸出信噪比?？梢钥闯?，QPT方法的輸出信噪比在相對速度為0時輸出信噪比與FrFT方法相同，隨著目標(biāo)相對速度提高，信號能量在變換域逐漸分散，輸出信噪比逐漸降低。上述數(shù)值分析表明，QPT方法的抗噪聲性能比FrFT方法差，該方法計算效率提高的代價是輸出信噪比的降低。雖然經(jīng)典的 FFT譜分析方法通過分析信號頻偏也可獲得目標(biāo)速度，但主動聲吶接收的目標(biāo)回波信號通常信噪比很低，F(xiàn)FT完全無法聚焦信號的能量，輸出信噪比遠(yuǎn)低于以上兩種方法。

圖11 不同相對速度QPT幅度譜的波形Fig.11 Waveforms of QPT spectrums corresponding to different relative velocities

圖12 QPT方法和FrFT方法的輸出信噪比對比Fig.12 Comparison of output signal to noise ratio between QPT method and FrFT method

4 結(jié) 論

本文提出了一種基于二次相位變換的主動聲吶線性調(diào)頻信號的速度估計方法，通過減少變換次數(shù)降低了 FrFT類算法速度估計的計算量。單次QPT速度估計方法利用LFM變換后幅度譜寬度與信號多普勒頻移的關(guān)系對速度進(jìn)行估計，速度符號的判斷依賴準(zhǔn)確的到達(dá)時刻。針對實際情況下難以準(zhǔn)確估計時延的問題，提出了對稱參數(shù)QPT速度符號判斷方法，通過對比對稱參數(shù)變換的幅度譜寬度判斷速度符號。由于不需要對變換參數(shù)進(jìn)行搜索，相較于 FrFT參數(shù)掃描的方法，所提速度估計方法可以降低一個維度的計算量。仿真實驗結(jié)果表明，該方法的速度估計精度可以達(dá)到1 m·s-1，如果只估計相對速度大小，計算量僅為精度相當(dāng)?shù)腇rFT參數(shù)掃描方法的1/60。但該方法計算效率提高的代價是輸出信噪比的降低，在實際應(yīng)用中需要進(jìn)行權(quán)衡?？紤]到復(fù)雜海洋環(huán)境條件下主動聲吶接收回波信號存在畸變，本文的方法還有待復(fù)雜信道下真實主動聲吶實驗數(shù)據(jù)的驗證。