丁 超 陳 喆 張宗堂 程玉勝
(海軍潛艇學院 青島 266000)
時延估計(Time Delay Estimation, TDE)是水聲領域的重要研究課題,基于時延估計的水聲目標被動測向、被動定位技術是水聲目標無源定位的重要分支。目前,時延估計通常可以用于以下3類算法:基于時延估計的目標方位估計算法[1],基于時延估計的目標被動定位算法[2],聯合其他觀測量(如到達角度、到達頻差等)的聯合定位方法[3,4]。在以上各類水聲目標被動測向和被動定位算法中,核心便是時延估計,時延估計的精度直接決定了測向和定位的精度。
本文結合互譜相位特征,設計了一種新的時延估計方法,在仿真實驗中取得了較好的處理結果。
建立信號模型為

其中,xi(t)為 接收器i所 接收到的信號,s (t)為目標輻射噪聲, τi為 s (t) 到達接收器i所 用的時間,αi(t)為s (t)到 達接收器i時 的衰減系數,ni(t)為 接收器i處的背景噪聲。理想假設下,n1(t), n2(t)與 s (t)互不相關。
將衰減系數αi(t)簡化為常量,為方便推導,在不影響模型本質的基礎上,將信號模型簡化為

其中, α為目標輻射噪聲到達兩接收器時的衰減系數比值, τ0為目標輻射噪聲到達兩接收器的時間差。所謂時延估計即是通過x1(t)和 x2(t)來 推算τ0的計 算過程。
目前,較為常用的時延估計方法有基本互相關法[5](Normalized Cross Correlation, NCC)、廣義互相關法[6,7](Generalized Cross Correlation,GCC)和自適應時延估計[8,9](Least Mean Square,LMS)等。
2.2.1 基本互相關算法
求取兩段波形最近似時的時延,經典算法便是基本互相關運算。結合式(2)信號模型,算式為

2.2.2 廣義互相關算法
實際情況中,目標輻射噪聲與背景噪聲間通常不是嚴格的互不相關,且兩路信號時長有限,背景噪聲會影響相關運算結果,干擾時延估計精度。因此需要在相關運算之前對信號進行濾波,以達到抑制噪聲、增強信號,從而提高時延估計精度的目的。結合式(2)信號模型,算式為


除以上算法之外,還可以利用互譜相位信息進行時延估計。結合式(2)信號模型,互譜算法為

通過求解互譜相位斜率,即可換算得到時延τ0[11–14]。假設存在單個目標且τ0=0.05s,理想情況下, φ(f) 為斜率2 πτ0的標準鋸齒波,如圖1(a)所示;非理想情況下,φ (f)為 斜率2 πτ0的疊加干擾的鋸齒波,為兩路寬帶信號添加信噪比10 dB的高斯白噪聲,信號時長5 s,采樣率10 kHz,其互譜相位如圖1(b)所示。
假設存在兩個目標且兩目標信號強度相同,為兩路寬帶信號添加信噪比10 dB的高斯白噪聲,信號時長5 s,采樣率10 kHz,其互譜相位如圖2所示。
此時互譜相位譜可以視為多個鋸齒波疊加干擾,可以表示為


圖1 單目標互譜相位圖

圖2 非理想假設下雙目標互譜相位圖
結合圖1和圖2,在互譜相位斜率求解過程中存在以下問題需要解決:
(1)互譜相位以2 π為周期變化,當相位變化超過2π 時,反正切函數的多值性會導致相位內卷。
針對該問題,目前較為常用的方法有相位拼接[15]、差分方法[16]等。相位拼接打破相位2 π的周期,通過周期間拼接,構造一條斜率2 πτ0的斜線;差分方法可將線性增量轉變為常量,從而進行線性增量估計。但相位拼接需要在相位拼接之前先進行一次是否需要相位修正的預判斷;差分方法需要在差分之后剔除相位在 ?π 與π 之間突變所產生的峰值,再進行常量估計。兩種方法均較為復雜,且準確性難以保證。
(2)受背景噪聲干擾,互譜相位出現波動甚至是跳變的現象,會導致斜率估計誤差增大。
針對該問題,目前較為常用的方法有最小二乘、等權平均[17]等。最小二乘法對相位數據進行線性擬合,以求得斜率參數;等權平均采用多點取樣并平均,以求得一個無偏的斜率估計。這兩種方法若想提高估算精度、提升抗噪能力,需要增加輸入點數、擴大輸入范圍,但在計算互譜相位斜率時,輸入范圍往往受限。
(3)當存在多目標時,互譜相位可能存在兩種情況:其一,形狀雜亂,無法計算斜率;其二,偶有斜率,是由多目標相位斜率疊加而成。這就導致目前各種直接計算相位斜率的算法都將不再適用。
正是由于以上困難,導致基于互譜相位斜率的時延估計方法難以用于實際復雜情況,本文正是針對以上困難設計了一種間接求取互譜相位斜率的算法,適用于強干擾、多目標的情況。
根據傅里葉變換的線性特性,結合式(8),可得多目標互譜相位的傅里葉變換為

其中,Si(k), N (k)分 別為si(f), n (f)的傅里葉變換。
可見,多目標互譜相位在頻域上相互疊加,加之噪聲干擾,難以直接求取相位斜率。但在進行傅里葉變換后,各目標互譜相位將以線譜形式分離,線譜位置能夠反映相位斜率大小,也就能進一步換算成時延大小。進行傅里葉變換后,通常的做法是對其取絕對值,但當輸入信號是互譜相位時,絕對值結果僅能反映相位斜率大小,不能反映相位斜率正負,也就無法確定時延正負,這會造成時延估計模糊。
根據周期信號傅里葉變換公式

由式(12)可知,鋸齒波傅里葉變換虛部在 f0處的值為a,其正負對應鋸齒波斜率的正負,如圖3所示。
可見,斜率為正的鋸齒波,其傅里葉變換虛部基頻處值為正;斜率為負的鋸齒波,其傅里葉變換虛部基頻處值為負。在互譜相位的處理中,互譜相位傅里葉變換虛部的線譜處橫坐標值即為兩路信號時延估計絕對值,線譜處數值的正負可以反映時延的正負。

圖3 鋸齒波及其傅里葉變換虛部
至此,依據互譜相位進行多目標時延估計的條件已經完備。本文據此設計出一種新的利用兩路信號互譜相位傅里葉變換域進行時延估計的算法,以互譜相位傅式變換法(Fourier transform domain of Cross-power Spectrum Phase, FCSP)表示,FCSP算法流程如圖4所示。
FCSP算法具有以下特點:
(1)利用互譜相位斜率進行時延估計,不會產生時延估計模糊;
(2)解決了現有的互譜相位斜率求解算法的諸多不便,不再受制于相位模糊,且能夠實現多目標分辨;
(3)通過多個互譜相位周期的累計,該算法具備較強的抗噪能力。

圖4 FCSP算法流程圖
LMS方法相當于一種迭代實現的GCC時延估計方法,效果與GCC方法一致,故本文僅對NCC,GCC與FCSP 3種算法展開比較,其中GCC算法采用常用的相位變換(PHAse Transform, PHAT)加權,即 ψ(f)=1/|Gx1,x2(f)|,該加權函數能夠銳化廣義互相關函數,突出時延峰值[18]。下面進行幾種情況下的仿真效果對比。
假設存在兩個目標且時延分別為0.05 s和0.03 s,為兩路寬帶信號添加信噪比–5 dB的高斯白噪聲,信號時長5 s,采樣率10 kHz,3種方法的時延估計結果如圖5所示。
為驗證算法性能,在3種算法計算結果歸一化之后,分別求取其噪聲的均方根誤差(Root Mean Squared Error, RMSE), RMSE越小則表示算法性能越好。進行3000次Monte Carlo仿真實驗,結果如表1所示。
由表1可知,FCSP算法結果線譜信噪比最高,最利于時延估計,GCC算法次之,NCC算法較差。
假設存在兩個目標且時延分別為0.05 s和0.03 s,目標信號包含50 Hz的強窄帶信號,為兩路寬帶信號添加信噪比–5 dB的高斯白噪聲,信號時長5 s,采樣率10 kHz,3種方法的時延估計結果如圖6所示。
由圖6可見,NCC算法無法消除強窄帶信號的干擾,當陣列最大時延大于窄帶信號周期時,便會出現時延估計模糊,對時延估計結果造成嚴重干擾。而GCC算法和FCSP算法可以消除時延估計模糊,且該情況下FCSP算法結果線譜信噪比更高。

圖5 多目標情況下3種算法性能比較

表1 3種算法噪聲RMSE均值對比

圖6 存在強窄帶信號情況下3種算法性能比較
假設存在兩個目標且時延分別為0.05 s和0.03 s,存在背景噪聲起伏,為兩路寬帶信號添加信噪比–5 dB的高斯白噪聲,信號時長5 s,采樣率10 kHz,3種方法的時延估計結果如7所示。
由圖7可見,NCC算法無法消除背景噪聲起伏的干擾,而GCC算法和FCSP算法可以消除,且該情況下FCSP算法結果線譜信噪比更高。
綜合3種情況來看,本文設計的FCSP算法能夠分辨多個目標,可以消除強窄帶信號及背景噪聲起伏的干擾,且具備較強的抗噪能力,在信噪比較高的情況下,輸出結果相較于GCC算法譜線信噪比更高。整體來看,FCSP算法處理效果優于NCC算法及GCC算法。
不同于目前常用的時延估計方法,本文從信號互譜相位角度出發設計了FGSP算法,利用兩路信號互譜相位的傅里葉變換域進行時延估計,并在仿真實驗中獲得了較好的處理結果。
(1)FGSP算法的設計是基于傅里葉變換的性質及互譜相位的特點,解決了一部分現有互譜相位斜率求解算法未解決的問題,準確求取互譜相位也就意味著準確求取了信號時延;

圖7 背景噪聲起伏情況下3種算法性能比較
(2)FGSP算法具備了多目標分辨能力、抗時延模糊能力、抗背景噪聲起伏能力及良好的抗噪能力,仿真效果整體優于NCC算法及GCC算法;
(3)FGSP算法的流程及思路不同于目前常用的NCC, GCC及LMS等算法,在后續的工作中,FGSP算法流程的各步驟還可進一步優化或者融入新方法,從而進一步提高時延估計精度。