一種互譜相位在時(shí)延估計(jì)中的應(yīng)用方法

丁 超 陳 喆 張宗堂 程玉勝

(海軍潛艇學(xué)院 青島 266000)

1 引言

時(shí)延估計(jì)(Time Delay Estimation, TDE)是水聲領(lǐng)域的重要研究課題，基于時(shí)延估計(jì)的水聲目標(biāo)被動(dòng)測向、被動(dòng)定位技術(shù)是水聲目標(biāo)無源定位的重要分支。目前，時(shí)延估計(jì)通?？梢杂糜谝韵?類算法：基于時(shí)延估計(jì)的目標(biāo)方位估計(jì)算法[1]，基于時(shí)延估計(jì)的目標(biāo)被動(dòng)定位算法[2]，聯(lián)合其他觀測量(如到達(dá)角度、到達(dá)頻差等)的聯(lián)合定位方法[3,4]。在以上各類水聲目標(biāo)被動(dòng)測向和被動(dòng)定位算法中，核心便是時(shí)延估計(jì)，時(shí)延估計(jì)的精度直接決定了測向和定位的精度。

本文結(jié)合互譜相位特征，設(shè)計(jì)了一種新的時(shí)延估計(jì)方法，在仿真實(shí)驗(yàn)中取得了較好的處理結(jié)果。

2 時(shí)延估計(jì)方法原理

2.1 時(shí)延估計(jì)基本模型

建立信號模型為

其中，xi(t)為 接收器i所 接收到的信號，s (t)為目標(biāo)輻射噪聲， τi為 s (t) 到達(dá)接收器i所 用的時(shí)間，αi(t)為s (t)到 達(dá)接收器i時(shí) 的衰減系數(shù)，ni(t)為 接收器i處的背景噪聲。理想假設(shè)下，n1(t), n2(t)與 s (t)互不相關(guān)。

將衰減系數(shù)αi(t)簡化為常量，為方便推導(dǎo)，在不影響模型本質(zhì)的基礎(chǔ)上，將信號模型簡化為

其中， α為目標(biāo)輻射噪聲到達(dá)兩接收器時(shí)的衰減系數(shù)比值， τ0為目標(biāo)輻射噪聲到達(dá)兩接收器的時(shí)間差。所謂時(shí)延估計(jì)即是通過x1(t)和 x2(t)來 推算τ0的計(jì) 算過程。

2.2 時(shí)延估計(jì)常用方法

目前，較為常用的時(shí)延估計(jì)方法有基本互相關(guān)法[5](Normalized Cross Correlation, NCC)、廣義互相關(guān)法[6,7](Generalized Cross Correlation,GCC)和自適應(yīng)時(shí)延估計(jì)[8,9](Least Mean Square,LMS)等。

2.2.1 基本互相關(guān)算法

求取兩段波形最近似時(shí)的時(shí)延，經(jīng)典算法便是基本互相關(guān)運(yùn)算。結(jié)合式(2)信號模型，算式為

2.2.2 廣義互相關(guān)算法

實(shí)際情況中，目標(biāo)輻射噪聲與背景噪聲間通常不是嚴(yán)格的互不相關(guān)，且兩路信號時(shí)長有限，背景噪聲會(huì)影響相關(guān)運(yùn)算結(jié)果，干擾時(shí)延估計(jì)精度。因此需要在相關(guān)運(yùn)算之前對信號進(jìn)行濾波，以達(dá)到抑制噪聲、增強(qiáng)信號，從而提高時(shí)延估計(jì)精度的目的。結(jié)合式(2)信號模型，算式為

3 互譜相位時(shí)延估計(jì)算法

3.1 互譜相位時(shí)延估計(jì)原理

除以上算法之外，還可以利用互譜相位信息進(jìn)行時(shí)延估計(jì)。結(jié)合式(2)信號模型，互譜算法為

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

假設(shè)存在兩個(gè)目標(biāo)且兩目標(biāo)信號強(qiáng)度相同，為兩路寬帶信號添加信噪比10 dB的高斯白噪聲，信號時(shí)長5 s，采樣率10 kHz，其互譜相位如圖2所示。

此時(shí)互譜相位譜可以視為多個(gè)鋸齒波疊加干擾，可以表示為

圖1 單目標(biāo)互譜相位圖

圖2 非理想假設(shè)下雙目標(biāo)互譜相位圖

3.2 現(xiàn)有互譜相位斜率求解方法及其限制

結(jié)合圖1和圖2，在互譜相位斜率求解過程中存在以下問題需要解決：

(1)互譜相位以2 π為周期變化，當(dāng)相位變化超過2π 時(shí)，反正切函數(shù)的多值性會(huì)導(dǎo)致相位內(nèi)卷。

針對該問題，目前較為常用的方法有相位拼接[15]、差分方法[16]等。相位拼接打破相位2 π的周期，通過周期間拼接，構(gòu)造一條斜率2 πτ0的斜線；差分方法可將線性增量轉(zhuǎn)變?yōu)槌Ａ?，從而進(jìn)行線性增量估計(jì)。但相位拼接需要在相位拼接之前先進(jìn)行一次是否需要相位修正的預(yù)判斷；差分方法需要在差分之后剔除相位在 ?π 與π 之間突變所產(chǎn)生的峰值，再進(jìn)行常量估計(jì)。兩種方法均較為復(fù)雜，且準(zhǔn)確性難以保證。

(2)受背景噪聲干擾，互譜相位出現(xiàn)波動(dòng)甚至是跳變的現(xiàn)象，會(huì)導(dǎo)致斜率估計(jì)誤差增大。

針對該問題，目前較為常用的方法有最小二乘、等權(quán)平均[17]等。最小二乘法對相位數(shù)據(jù)進(jìn)行線性擬合，以求得斜率參數(shù)；等權(quán)平均采用多點(diǎn)取樣并平均，以求得一個(gè)無偏的斜率估計(jì)。這兩種方法若想提高估算精度、提升抗噪能力，需要增加輸入點(diǎn)數(shù)、擴(kuò)大輸入范圍，但在計(jì)算互譜相位斜率時(shí)，輸入范圍往往受限。

(3)當(dāng)存在多目標(biāo)時(shí)，互譜相位可能存在兩種情況：其一，形狀雜亂，無法計(jì)算斜率；其二，偶有斜率，是由多目標(biāo)相位斜率疊加而成。這就導(dǎo)致目前各種直接計(jì)算相位斜率的算法都將不再適用。

正是由于以上困難，導(dǎo)致基于互譜相位斜率的時(shí)延估計(jì)方法難以用于實(shí)際復(fù)雜情況，本文正是針對以上困難設(shè)計(jì)了一種間接求取互譜相位斜率的算法，適用于強(qiáng)干擾、多目標(biāo)的情況。

3.3 互譜相位斜率求解方法

根據(jù)傅里葉變換的線性特性，結(jié)合式(8)，可得多目標(biāo)互譜相位的傅里葉變換為

其中，Si(k), N (k)分 別為si(f), n (f)的傅里葉變換。

可見，多目標(biāo)互譜相位在頻域上相互疊加，加之噪聲干擾，難以直接求取相位斜率。但在進(jìn)行傅里葉變換后，各目標(biāo)互譜相位將以線譜形式分離，線譜位置能夠反映相位斜率大小，也就能進(jìn)一步換算成時(shí)延大小。進(jìn)行傅里葉變換后，通常的做法是對其取絕對值，但當(dāng)輸入信號是互譜相位時(shí)，絕對值結(jié)果僅能反映相位斜率大小，不能反映相位斜率正負(fù)，也就無法確定時(shí)延正負(fù)，這會(huì)造成時(shí)延估計(jì)模糊。

根據(jù)周期信號傅里葉變換公式

由式(12)可知，鋸齒波傅里葉變換虛部在 f0處的值為a，其正負(fù)對應(yīng)鋸齒波斜率的正負(fù)，如圖3所示。

可見，斜率為正的鋸齒波，其傅里葉變換虛部基頻處值為正；斜率為負(fù)的鋸齒波，其傅里葉變換虛部基頻處值為負(fù)。在互譜相位的處理中，互譜相位傅里葉變換虛部的線譜處橫坐標(biāo)值即為兩路信號時(shí)延估計(jì)絕對值，線譜處數(shù)值的正負(fù)可以反映時(shí)延的正負(fù)。

圖3 鋸齒波及其傅里葉變換虛部

至此，依據(jù)互譜相位進(jìn)行多目標(biāo)時(shí)延估計(jì)的條件已經(jīng)完備。本文據(jù)此設(shè)計(jì)出一種新的利用兩路信號互譜相位傅里葉變換域進(jìn)行時(shí)延估計(jì)的算法，以互譜相位傅式變換法(Fourier transform domain of Cross-power Spectrum Phase, FCSP)表示，F(xiàn)CSP算法流程如圖4所示。

FCSP算法具有以下特點(diǎn)：

(1)利用互譜相位斜率進(jìn)行時(shí)延估計(jì)，不會(huì)產(chǎn)生時(shí)延估計(jì)模糊；

(2)解決了現(xiàn)有的互譜相位斜率求解算法的諸多不便，不再受制于相位模糊，且能夠?qū)崿F(xiàn)多目標(biāo)分辨；

(3)通過多個(gè)互譜相位周期的累計(jì)，該算法具備較強(qiáng)的抗噪能力。

4 算法效果仿真比較

圖4 FCSP算法流程圖

LMS方法相當(dāng)于一種迭代實(shí)現(xiàn)的GCC時(shí)延估計(jì)方法，效果與GCC方法一致，故本文僅對NCC,GCC與FCSP 3種算法展開比較，其中GCC算法采用常用的相位變換(PHAse Transform, PHAT)加權(quán)，即 ψ(f)=1/|Gx1,x2(f)|，該加權(quán)函數(shù)能夠銳化廣義互相關(guān)函數(shù)，突出時(shí)延峰值[18]。下面進(jìn)行幾種情況下的仿真效果對比。

4.1 多目標(biāo)情況

假設(shè)存在兩個(gè)目標(biāo)且時(shí)延分別為0.05 s和0.03 s，為兩路寬帶信號添加信噪比–5 dB的高斯白噪聲，信號時(shí)長5 s，采樣率10 kHz，3種方法的時(shí)延估計(jì)結(jié)果如圖5所示。

為驗(yàn)證算法性能，在3種算法計(jì)算結(jié)果歸一化之后，分別求取其噪聲的均方根誤差(Root Mean Squared Error, RMSE), RMSE越小則表示算法性能越好。進(jìn)行3000次Monte Carlo仿真實(shí)驗(yàn)，結(jié)果如表1所示。

由表1可知，F(xiàn)CSP算法結(jié)果線譜信噪比最高，最利于時(shí)延估計(jì)，GCC算法次之，NCC算法較差。

4.2 存在強(qiáng)窄帶信號情況

假設(shè)存在兩個(gè)目標(biāo)且時(shí)延分別為0.05 s和0.03 s，目標(biāo)信號包含50 Hz的強(qiáng)窄帶信號，為兩路寬帶信號添加信噪比–5 dB的高斯白噪聲，信號時(shí)長5 s，采樣率10 kHz，3種方法的時(shí)延估計(jì)結(jié)果如圖6所示。

由圖6可見，NCC算法無法消除強(qiáng)窄帶信號的干擾，當(dāng)陣列最大時(shí)延大于窄帶信號周期時(shí)，便會(huì)出現(xiàn)時(shí)延估計(jì)模糊，對時(shí)延估計(jì)結(jié)果造成嚴(yán)重干擾。而GCC算法和FCSP算法可以消除時(shí)延估計(jì)模糊，且該情況下FCSP算法結(jié)果線譜信噪比更高。

圖5 多目標(biāo)情況下3種算法性能比較

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

圖6 存在強(qiáng)窄帶信號情況下3種算法性能比較

4.3 背景噪聲起伏情況

假設(shè)存在兩個(gè)目標(biāo)且時(shí)延分別為0.05 s和0.03 s，存在背景噪聲起伏，為兩路寬帶信號添加信噪比–5 dB的高斯白噪聲，信號時(shí)長5 s，采樣率10 kHz,3種方法的時(shí)延估計(jì)結(jié)果如7所示。

由圖7可見，NCC算法無法消除背景噪聲起伏的干擾，而GCC算法和FCSP算法可以消除，且該情況下FCSP算法結(jié)果線譜信噪比更高。

綜合3種情況來看，本文設(shè)計(jì)的FCSP算法能夠分辨多個(gè)目標(biāo)，可以消除強(qiáng)窄帶信號及背景噪聲起伏的干擾，且具備較強(qiáng)的抗噪能力，在信噪比較高的情況下，輸出結(jié)果相較于GCC算法譜線信噪比更高。整體來看，F(xiàn)CSP算法處理效果優(yōu)于NCC算法及GCC算法。

5 結(jié)束語

不同于目前常用的時(shí)延估計(jì)方法，本文從信號互譜相位角度出發(fā)設(shè)計(jì)了FGSP算法，利用兩路信號互譜相位的傅里葉變換域進(jìn)行時(shí)延估計(jì)，并在仿真實(shí)驗(yàn)中獲得了較好的處理結(jié)果。

(1)FGSP算法的設(shè)計(jì)是基于傅里葉變換的性質(zhì)及互譜相位的特點(diǎn)，解決了一部分現(xiàn)有互譜相位斜率求解算法未解決的問題，準(zhǔn)確求取互譜相位也就意味著準(zhǔn)確求取了信號時(shí)延；

圖7 背景噪聲起伏情況下3種算法性能比較

(2)FGSP算法具備了多目標(biāo)分辨能力、抗時(shí)延模糊能力、抗背景噪聲起伏能力及良好的抗噪能力，仿真效果整體優(yōu)于NCC算法及GCC算法；

(3)FGSP算法的流程及思路不同于目前常用的NCC, GCC及LMS等算法，在后續(xù)的工作中，F(xiàn)GSP算法流程的各步驟還可進(jìn)一步優(yōu)化或者融入新方法，從而進(jìn)一步提高時(shí)延估計(jì)精度。