基于頻譜Sevcik維數的目標運動速度估計*

趙曉靜 高大治 孫 凱 李玉征 郭豆豆

(中國海洋大學信息科學與工程學部 青島 266100)

0 引言

線譜具有能量高、傳播距離遠等特點[1]，被廣泛應用于目標運動參數估計。線譜多普勒頻移包含著目標的運動速度、到達最近點時間、最近距離等運動參數信息。1994年，Ferguson 等[2]提出了利用線譜多普勒頻移估計目標運動參數的方法，該方法利用時頻分析獲取瞬時頻率，與多普勒頻移模型進行擬合，獲取目標的運動參數。而后，Quinn[3]推導了接收器和運動聲源運動方向不共線時多普勒頻移表達式。鄒紅星等[4]提出了Dopplerlet 變換，該變換可將多普勒頻移去除。吳國清等[5]利用鄒紅星等提出的Dopplerlet 變換和匹配投影算法，提取線譜多普勒中包含的目標速度和最近距離等運動信息，而后又提出了一種利用多線譜的多普勒頻移估計運動目標最近距離的方法[6]，先用Wigner-Ville分布提取瞬時頻率，再定義多普勒頻移的基函數，實現目標的運動參數估計。該方法的缺點是需要大量的先驗知識。高偉等[7]提出了一種基于目標噪聲輻射強度和線譜多普勒頻移估計目標最近距離的方法，該方法不需要大量的先驗知識，且搜索范圍易確定。徐靈基[8]針對水下目標線譜多普勒進行深入研究，通過引入了一種時頻變換–匹配Wigner變換，獲取多普勒頻移，從而實現目標運動參數估計。Liang[9]等提出了Doppler-Chirplet 變換，該變換直接利用觀測到的瞬時頻率與源速度的關系作為Chirplet 變換的核心，可直接估測運動目標的速度。同時，還有多尺度Chirplet 路徑追蹤方法[10]等。劉凱悅等[11]提出了一種基于線譜多普勒的水下對空中聲源的運動參數獲取方法。這些算法不僅提高了計算精度，也縮短了計算時間。現有的基于多普勒頻移獲取運動參數大都基于瞬時頻率。高德洋等[12]提出了一種不需要提取瞬時頻率的線譜目標運動參數獲取方法，他首先提出Doppler-warping變換，將warping 變換的思想運用到多普勒頻移的去除上，又利用頻譜求熵的方法獲取目標的運動參數，相較于Ferguson 等的算法，不需要提取瞬時頻率，因此適用性更高，計算精度也更高。

分形維數與熵一樣，都是表征系統復雜程度的物理量。因此，分形維數在信號處理方面有廣泛的應用?，F階段，分形維數被廣泛應用于圖像分類[13-14]、故障診斷[15]、地形分析[16]、頻譜感知[17-18]等方面。截至目前，鮮有學者將分形維數應用于線譜信號目標運動參數估計中。本文將頻譜分形維數應用到目標運動速度獲取中。首先利用Doppler-warping變換對接收信號重采樣，然后通過計算頻譜的分形維數作為衡量多普勒頻移消除程度的代價函數，最后通過獲取代價函數最小值對應的速度值獲取目標的運動速度。

1 理論

1.1 多普勒頻移及Doppler-warping變換

當聲源沿不經過接收器的路徑做勻速直線運動時，接收器接收信號存在多普勒頻移，其運動示意圖如圖1 所示。其中v為目標聲源的運動速度，t0為到達最近點時間，r0為最近距離，接收信號瞬時頻率由Quinn[3]推導的表達式得出，如式(1) 所示：

圖1 聲源與接收器運動示意圖Fig.1 Geometrical model of the wayside acoustic testing model

其中，f0為目標聲源發射信號的頻率。本文中，設置t0、r0已知，搜索目標勻速直線運動速度v。由式(1)可以看出，由于多普勒效應的存在，接收信號的頻率不再為單頻，而是隨時間變化的。

Doppler-warping 變換由中國海洋大學水聲研究團隊[12]提出，warping 算子的表達式如式(2)所示：

使用式(2)對接收信號進行重采樣，即將式(2)代入式(1)的t中，經整理后f(t)=f0，具體推導過程見文獻[12]。此時多普勒頻移被消除，原始信號的單頻特性恢復。接收信號經過Doppler-warping 變換前后，LOFAR圖如圖2所示。

圖2 接收信號Doppler-warping 變換前后Fig.2 LOFAR figure of single-frequency signal and LOFAR figure after Doppler-warping transform

1.2 Sevcik維數

分形維數是衡量事物分形特性的一種度量參數，也是衡量事物復雜程度的物理量。分形維數主要可分為Katz 分形維數、Hausdorff維數、計盒維數、Sevcik 分形維數和自相似分形維數等幾種。其中，Sevcik 維數多用于衡量信號的波形的分形特性，且計算過程較簡單，計算時間短。因此，本文采用Sevcik 維數來衡量信號頻譜的分形特性。Sevcik 維數的表達式如式(3)所示：

其中，N為信號點個數，L為信號波形長度，其表達式如式(4)所示：

對于一個時域信號s(t)，若其時長為T，采樣率為fs，信號點個數為N，則其經過快速傅里葉變換(Fast Fourier transform,FFT)后，頻譜S(f)的xk對應信號頻率點fk。此時其信號波形長度L的表達式經整理后如式(6)所示：

1.3 算法步驟

由于多普勒的存在，線譜信號的頻譜能量不再集中在原始頻率上，而是在中心頻率周圍有一定的展寬，此時，信號頻譜的分形維數會比不存在多普勒時明顯升高，通過Doppler-warping 變換后，可將多普勒頻移去除，此時信號恢復單頻特性，其頻譜的分形維數會明顯降低。Doppler-warping 算子中包含最近點時間、最近距離、速度等運動參數，因此，Doppler-warping算子中運動參數的準確性影響多普勒的去除程度?？赏ㄟ^線譜和連續譜結合的方式，估計目標到達最近點的時間和最近距離，本文假設以上兩參數已知，利用Doppler-warping 變換結合頻譜Sevcik 維數單獨進行速度估計，其步驟可概括為以下幾步：

(1)設置速度搜索網格[v1,v2,···,vk]。應用時，可根據需要設置速度搜索步長，在常見目標速度估計應用中可將搜索步長設為0.1 m/s。

(2) 最近點時間和最近距離兩參數由其他方法估計，設為已知，計算每個速度值下的Dopplerwarping算子。

(3) 將不同算子代入接收時域信號s(t)，對接收信號進行重采樣，獲得信號s′(t)。

(4) 對s′(t)進行傅里葉變換后獲得頻譜S′(f)，取模并計算每一個S′(f)的Sevcik分形維數D(vk)。

(5) 找到D(vk)取值最小時對應的速度v，即為速度估計值。

算法流程圖如圖3所示。

圖3 算法流程Fig.3 Flow chart of the algorithm

2 仿真及算法性能

2.1 仿真及結果

基于以上理論，進行仿真驗證，本文算法對空氣中目標和水中目標同樣適用，這里以空氣中目標為例。聲源信號為600 Hz 單頻信號，采樣率為10000 Hz，信號時長10 s。聲源運動速度為10 m/s，最近距離為10 m，最近點時間為5 s。聲速為340 m/s。接收信號及處理結果如圖4所示。

圖4 仿真結果Fig.4 Simulation results

設置速度搜索網格為[1 : 0.2 : 20]，搜索速度v，結果如圖4(d)所示。

利用本文方法進行搜索速度，從圖4(d)可以看出，在真實速度值10 m/s 處，Sevcik 維數明顯低于錯誤值處，因此能準確搜索出速度結果，證明本文提出的方法有效。

2.2 算法性能

為了體現本文提出算法的性能，與用熵作為代價函數的搜索結果進行比較。在相同的仿真條件下，分別用熵作為代價函數和Sevcik 維數作為代價函數，進行速度估計。估計結果如圖5所示。

圖5 熵和Sevcik 搜索結果對比Fig.5 Comparision of estimation result of entropy and Sevcik

由于本文算法是尋找代價函數最小值，因此定義最小代價函數值為0 dB，定義3 dB 處對應的兩速度值的差為主瓣寬度。根據以上定義，用Sevcik維數作為代價函數時，主瓣寬度為0.76 m/s，而用熵作為代價函數，主瓣寬度為3.81 m/s。可見用Sevcik維數和熵作為代價函數時，雖然都可以搜索出正確速度值，但是本文提出的算法，搜索結果在正確值位置處起伏更明顯，瓣更窄，性能更好。

在實際應用中，往往接收信號存在噪聲，有時甚至信噪比(Signal to noise ratio,SNR) 很低，因此，本文利用仿真分析不同SNR下兩種算法的性能對比。仿真條件為：信號時長20 s，聲源運動速度為10 m/s，最近距離為10 m，最近點時間為10 s。分別在仿真的接收信號中加入高斯白噪聲，SNR計算方法采用譜SNR。兩種算法搜索結果對比如圖6 和表1所示。

表1 典型SNR 下兩方法速度估計誤差對比Table 1 Comparison of speed estimation errors between the two methods under typical SNR

圖6 SNR 對兩種算法搜索結果的影響Fig.6 Influence of SNR on the two algorithms

從圖6 和表1 中可以看出，在仿真條件相同的情況下，SNR 過低時，Sevcik 維數和熵作為代價函數效果均不佳，這是因為當SNR 較低時，噪聲對信號頻譜有較大影響，因此兩種方法均會失效。但是，從圖6 可以看出，Sevcik 維數作為代價函數抗噪性提高5 dB 左右。主要原因是兩種算法都基于頻譜的混亂程度，即使在SNR 低于-5 dB 的情況下，多普勒的存在與否也會對頻譜的Sevcik 維數產生較大影響，即存在多普勒時頻譜的Sevcik 維數明顯大于不存在時，因此本文算法可行。而只有在SNR大于-5 dB 以后，用熵值作為代價函數時，多普勒的存在與否對頻譜的影響較大，熵值作為代價函數的算法才會有效。SNR 對多普勒是否存在時兩種不同代價函數的影響如圖7(a)、圖7(b)所示。在相同SNR下，多普勒的存在對頻譜Sevcik維數的影響明顯大于頻譜熵(如圖7(c)所示)。

圖7 SNR 對多普勒是否存在時兩種代價函數的影響Fig.7 The effect of the SNR on the presence or absence of Doppler on the two cost functions

3 實驗及結果

3.1 實驗

針對上述理論和仿真結果，進行實驗。陸面目標采用某卡車運動過程中輻射的線譜噪聲，數據來源為acousticstoday.org。該信號時長30 s，采樣頻率為12000 Hz。聲速為347 m/s。最近距離為35 m，速度為20 km/h。其到達最近點時間為15.2 s，本文在處理過程中將其設為已知量。接收信號LOFAR圖以及搜索結果如圖8所示。

圖8 卡車實驗結果Fig.8 Truck experiment results

圖8(a)為接收器接收信號LOFAR 圖，該圖中線譜多普勒十分明顯。設置速度網格為[1:0.1:15]，單獨搜索速度的結果如圖8(b)所示，結果為5.8 m/s，相對誤差僅為4.4%，該結果證明了本文算法的有效性。本文方法與熵值方法搜索結果對比如圖8(c)所示。計算兩種算法速度搜索結果的主瓣寬度，熵值算法的主瓣寬度為3.15 m/s，本文算法的主瓣寬度僅為0.74 m/s，結果表明本文算法主瓣寬度明顯更窄。

為了進一步證明本文算法的有效性，選取2022年6 月在青島近海進行的海試實驗數據，實驗聲源為航道上經過的貨輪，該目標為非合作目標，貨輪的運動軌跡經接收船船載AIS(船舶自動識別系統)記錄如圖9(a)所示，距離變換如圖9(b)所示。接收信號為貨輪在運動過程中輻射線譜噪聲。實驗地聲速為1510 m/s，選取的信號時長為1200 s，信號采樣頻率為80000 Hz。由于信號采樣率過高，而處理信號所選線譜頻段僅為980～1000 Hz 左右，因此對接收信號進行降采樣以減少計算量。接收信號頻譜圖和LOFAR 圖如圖9(c)、圖9(d) 所示。貨輪運動平均速度為5.3 m/s，到達水聽器最近距離為2010 m，選取的信號中，目標到達最近點的時間為49 s。利用本文提出方法進行運動速度搜索，搜索結果如圖9(e)所示。設置速度網格為[1:0.1:10]，速度搜索結果為5.1 m/s，相對誤差僅為3.8%。本文算法和熵值算法速度搜索對比結果如圖9(f)所示。計算兩種算法的主瓣寬度，熵值算法的主瓣寬度為1.31 m/s，而本文算法的主瓣寬度僅為0.28 m/s，因此本文算法的性能較熵值算法有明顯提高。

圖9 海試實驗結果Fig.9 Marine experiment results

以上兩組實驗結果誤差分別為4.4%和3.8%，均能較為準確估計目標運動速度值，速度搜索相對誤差均小于5%。兩組實驗結果均證明了本文提出算法的有效性。

4 結論

本文提出了Doppler-warping 變換結合頻譜分形維數的方法來獲取目標的運動速度。傳統方法通過迭代算法提取接收信號的瞬時頻率，計算時間較長，本文提出的算法過程簡單，計算時間較短。相對于用熵作為衡量頻譜復雜度的代價函數，本文利用Sevcik維數作為代價函數，經仿真檢驗，抗噪能力較熵值方法提高5 dB。經仿真和實驗檢驗，本文算法代價函數主瓣寬度更窄，具有更好的性能?？諝夂退新暷繕藢嶒灳蓽蚀_搜索出目標運動速度，兩實驗誤差均小于5%，證明了本文所提算法的有效性。下一步擬將本文方法用于寬帶噪聲目標和雙目標的運動參數獲取中，實現算法應用的擴展。