基于水聲擴頻信號的空時分簇DoA 估計算法

周鋒，張寶勝，張文博

（1.哈爾濱工程大學水聲技術重點實驗室，黑龍江 哈爾濱 150001；2.海洋信息獲取與安全工信部重點實驗室（哈爾濱工程大學），黑龍江 哈爾濱 150001；3.哈爾濱工程大學水聲工程學院，黑龍江 哈爾濱 150001）

0 引言

在復雜的海洋環(huán)境中，水聲擴頻通信具有抗多徑、抗干擾、多址通信等優(yōu)點，是實現(xiàn)高質量水聲通信的首選技術之一[1]。為了進一步提高系統(tǒng)性能，研究人員提出了利用水聽器陣列接收信號。將陣列信號處理技術應用于水聲擴頻通信，首先要解決的是多徑信號波達方向（DoA,direction of arrival）估計問題，因為DoA 信息是完成波束成形、時空二維Rake 接收的重要依據(jù)，所以本文主要研究水聲擴頻信號的DoA 問題。

水聲擴頻信號是寬帶信號，目前比較經典的寬帶高分辨算法主要有2 類[2]，分別是非相干信號子空間方法（ISM,incoherent signal subspace method）[3]和相干信號子空間方法（CSM,coherent signal subspace method）[4]。前者利用頻率分解的思想，將寬帶源信號分解為一系列窄帶信號再進行DoA 估計。后者基于頻率聚焦的原理，其中最典型的為雙邊相關變換（TCT,two-sided correlation transformation）方法[5]，利用聚焦矩陣將分解的子信號變換到參考頻率上，進而使用窄帶處理方法實現(xiàn)DoA 估計。但這樣的處理方法使CSM 需要預估信號源角度，而聚焦矩陣對預估角度的依賴使最終DoA 估計結果產生偏差。

水下通信信號的DoA 估計往往需要考慮比空氣中通信信號的DoA 估計更多的因素，水下的地形崎嶇，聲波在水中傳播容易出現(xiàn)散射、折射，造成陣列接收信號包含相干信號，可能導致DoA 估計算法失效[6]。由于海水流動或者通信雙方的相對運動，信道結構發(fā)生變化，導致接收信號的DoA和信源個數(shù)隨之變化。同時，擴頻通信常常工作在較低信噪比下，使許多常規(guī)的DoA 估計算法失效。

針對上述問題，研究人員在水下DoA 估計領域展開了大量研究。劉亮等[7]提出了一種適用于水聲寬帶信號的DoA 估計算法，對每個頻點信號進行不同的加權，在低信噪比下獲得了更好的分辨率。韓澤洋等[8]針對多徑問題，提出了一種線性調頻信號波達方向估計方法。文獻[9]提出了一種基于進化算法思想的水下DoA 估計算法。文獻[10]提出了一種基于信號自消零的低快拍數(shù)水下DoA 估計算法。文獻[11]利用m 序列的天然抗多徑的能力，使各條路徑相干性減弱，達到天然解相干的效果，解決了多重信號分類（MUSIC,multiple signal classification）算法中相干源的問題。

然而，以上方法并沒有考慮DoA 估計的抗干擾問題，將干擾信號的DoA 信息作為波束導向會導致水聲擴頻信號被削弱，對通信系統(tǒng)造成負面影響。為此，本文提出一種基于水聲擴頻信號的空時分簇（STDC,space-tim division cluster）DoA 估計算法，該算法不需要預先判斷信源個數(shù)，不受干擾影響，只估計水聲擴頻信號的DoA，且在低信噪比下依然保持較高的角度分辨率，具有較好的穩(wěn)健性。

1 信號模型

假設發(fā)射信號是一段由m 序列對單頻信號進行擴頻的信號，表示為

本文以均勻直線陣為例，假設聲源在遠場，經過信道后，有P條路徑上的水聲擴頻信號入射到M元均勻直線陣列上，圖1 為第p條路徑以角度θp入射到陣列上。

圖1 水聲擴頻信號入射到均勻直線陣列

其中，第一個陣元為參考陣元，d為陣元間距，則第m個陣元接收的信號為

其中，a為多普勒壓縮因子，τp為第p條路徑到達參考陣元的時延，τm(θp)為第m個陣元相對于參考陣元的聲程差，nm(t)為第m個陣元的噪聲。

2 空時分簇DoA 算法原理與分析

為了穩(wěn)健地估計水聲擴頻信號的入射角度，本文設計了一種基于水聲擴頻信號的空時分簇DoA估計算法，算法框架如圖2 所示。首先，對各個陣元接收的信號進行正交解調，對解調后的數(shù)據(jù)進行多波束成形，生成 ?90°～90°的波束成形后的數(shù)據(jù)；然后，對波束成形后的數(shù)據(jù)進行時頻二維搜索，消除多普勒效應的影響；最后，搜索出最大值，得到信號的入射角度θ、頻偏Δω和時延τ。

圖2 基于水聲擴頻信號的空時分簇DoA 估計算法框架

2.1 正交解調

由于信號在基帶上便于處理，接收到的信號通常下變頻到基帶。如圖3 所示，將第m個陣元接收的信號乘以余弦和正弦信號，并通過低通濾波器（LPF,low pass filter），這個過程稱為正交解調[12]。

圖3 正交解調

正交解調的輸出是由同向分量和正交分量組成的復數(shù)，即

其中，nmIQ(t)=nmI(t)+jnmQ(t)，且j=。

2.2 多波束成形

目前，有許多波束成形方法，為了簡化方法，本文使用相移波束成形[13]。相移波束成形的原理是將每一個陣元的信號以相同的相位疊加，如圖4 所示。

圖4 相移波束成形原理

此時，每個陣元的相位差是由符號相位和載波相位決定的。為了實現(xiàn)相移波束成形，必須先對擴頻序列的符號相位進行對齊。第m個陣元與參考陣元的聲程差為

利用聲程差進行符號相位補償后，第m個陣元的信號為

設第m個陣元的權值為ωm，則將所有陣列元素的加權輸出相加，作為波束成形后的陣列輸出，即

其中

設置陣元個數(shù)為10，陣元間距為18 cm，聲速為1 500 m/s，載波頻率為4 000 Hz，上述波束成形的靜態(tài)方向圖如圖5 所示。

圖5 波束成形的靜態(tài)方向圖

當計算量不滿足要求時，可通過增加波束掃描間隔來減少計算量。上述波束成形方法的波束寬度與常規(guī)波束寬度類似，故借用常規(guī)波束寬度的計算式，圖1 所示的均勻直線陣列的波束寬度BW 為

其中，信號波長λ=，c為聲速，fc為載波中心頻率。

多波束成形的原理是采用不同組的相位補償獲得多個不同方向的波束，各波束間應有一定的覆蓋，否則會造成目標遺漏[14]。通常以波束寬度一半處為覆蓋點，多波束間的覆蓋如圖6 所示。

圖6 多波束間的覆蓋

2.3 時頻二維搜索

根據(jù)文獻[11]，對正交解調后的復基帶信號進行時頻二維搜索，可得到信號的時延和多普勒頻偏。時頻二維搜索如圖7 所示，首先利用不同的頻偏對基帶信號進行補償，然后利用本地的擴頻序列進行互相關，最后得到用于估計多普勒頻偏和時延的模糊函數(shù)，其中 Δω表示頻偏。

圖7 時頻二維搜索

如果本地擴頻序列很長，那么完成一次相關需要很長時間，不能保證通信的實時性要求。但是如果考慮到時域卷積等于頻域相乘，那么可以將時域相關運算轉化到頻域做相乘運算，這樣可以大大減少相關所需的時間[15]。其基本原理為

其中，y(n)表示經過頻偏補償和離散后的基帶信號，c(n)表示本地的擴頻序列，?表示卷積，F(xiàn)FT表示快速傅里葉變換，F(xiàn)FT?表示快速傅里葉變換的共軛，IFFT 表示快速傅里葉逆變換，F(xiàn)FT[y(n)]表示將時域信號y(n)變換到頻域，F(xiàn)FT*[c(n)]表示將時域信號c(n)變到頻域然后取共軛?；贔FT 的快速相關算法框架如圖8 所示，其中，|?|2表示模值的平方。

圖8 基于FFT 的快速相關算法框架

頻偏為10 Hz 時7 階m 序列的時頻二維搜索仿真如圖9 所示。

圖9 頻偏為10 Hz 時7 階m 序列的時頻二維搜索仿真

從圖9 可以看出，搜索峰具有一定的寬度，說明頻偏估計具有一定的分辨率，與波束掃描一樣，可適當增加搜索間隔減少計算量，且采用不同長度的擴頻序列其估計精度將不同。假設本地的擴頻序列與信號中的擴頻序列對齊且相乘時為1，忽略噪聲和多徑干擾的影響，對式(8)中的第p條路徑的波束成形后的水聲擴頻信號做相關運算并平方后可得

其中，Ts為采樣時間，N為擴頻序列的周期，A為信號幅度。至此，可以通過式(13)看到載波頻偏對水聲擴頻信號的相關峰值造成的影響，且擴頻序列的長度不同，式(13)的主瓣寬度也不同。水聲擴頻信號頻率最小分辨率 Δf為[16-17]

為了提高頻率搜索的速度，可以設置合適的搜索頻偏步長以減少計算量，7 階m 序列的頻率分辨率如圖10 所示。

圖10 7 階m 序列的頻率分辨率

尋找空時頻三維搜索中的最大值可以得到接收水聲擴頻信號的頻偏、時延和角度，按照估計所得頻偏，將補償后沒有多普勒頻偏的相關輸出挑選出來，便可以得到入射信號的空時分簇圖。

3 仿真驗證

本文以水聲擴頻通信為背景，估計水聲擴頻信號的DoA，目的在于使擴頻通信中DoA 估計算法在低信噪比以及復雜的環(huán)境下具有穩(wěn)定的性能。

在理想情況（無編碼擴頻）下，7 階m 序列的擴頻通信相比于QPSK 有21 dB 的SNR 增益，陣元個數(shù)為10 個，波束成形后相對于無編碼擴頻的SNR增益為10 dB，如圖11 所示。以下仿真都以此信噪比作為基準。

圖11 7 階m 序列的擴頻通信誤比特率對比

實驗條件為10 個間距d=0.18 m的陣元組成的均勻直線陣列，聲速為1500 m/s，水聲擴頻信號中心頻率為4 kHz，帶寬為2.5 kHz。

在MATLAB 仿真環(huán)境下分析本文所提空時分簇DoA 估計算法的性能，對比算法為傳統(tǒng)波束成形（CBF,conventional beamforming）算法、文獻[11]中的MUSIC 算法、ISM、CSM-TCT。假設以上子空間類算法的信源數(shù)已知，ISM 采用前向空間平滑進行解相干，CSM-TCT 算法中預估角度已知，所提DoA 估計算法的使用條件是提前已知發(fā)射信號的擴頻序列，但不需要知道信源數(shù)和角度預估值。

3.1 抗干擾性分析

假設水聲擴頻信號的入射角為30°，有一個頻率為4 000 Hz 的單頻干擾信號的入射角度為?10°，信噪比為?10 dB，信干比為?10 dB，其他條件保持不變。空時分簇DoA 估計算法的仿真結果如圖12 所示，對比算法的仿真結果如圖13 所示。

圖12 入射角為30°時空時分簇DoA 估計算法的仿真結果

圖13 入射角為30°時對比算法的仿真結果

本文所提空時分簇DoA 估計算法利用本地的擴頻序列與多波束成形后的接收信號進行循環(huán)相關，可以充分利用擴頻序列的抗干擾性，抑制干擾，且只估計水聲擴頻信號的方向，其他干擾信號的方向均不顯示。

對比圖12 和圖13 可知，空時分簇DoA 估計算法不會受到干擾信號的影響，而MUSIC 算法、CBF 算法、ISM、CSM-TCT 在強干擾的影響下，干擾信號的主瓣比水聲擴頻信號的主瓣高，若選取最大值作為有用信號的入射方向，將會影響通信，而空時分簇DoA 估計可以充分發(fā)揮水聲擴頻信號的抗干擾優(yōu)勢。

3.2 分辨率分析

假設水聲擴頻信號經過信道后，有2 條主要路徑到達接收陣列，信號的入射角分別為30°和25°，信噪比為?10 dB?？諘r分簇DoA 估計算法的仿真結果如圖14 所示，對比算法的仿真結果如圖15 所示。

圖14 入射角分別為30°和25°時空時分簇DoA 估計算法的仿真結果

圖15 入射角分別為30°和25°時對比算法的仿真結果

信噪比為?10 dB 時，4 種對比算法已經無法分辨2 個入射角度，而所提空時分簇DoA 估計算法可以分辨2 個入射角度。這是因為信號經過信道后經過的路徑不同，所以傳播時間也會有差別，可以利用時間和角度將每條路徑分別開。

3.3 抗噪聲分析

本節(jié)對比不同算法的抗噪聲性能。假設水聲擴頻信號以30°入射到天線陣列，信噪比SNR 為?35～10 dB，不同算法的估計均方根誤差隨信噪比變化如圖16 所示。

圖16 不同算法的估計均方根誤差隨信噪比變化

從圖16 可以看出，隨著信噪比的增加，幾種算法的均方根誤差都隨之下降。相同信噪比下，空時分簇DoA 估計算法的性能最好。在工程誤差允許的范圍內，空時二維搜索的最低工作信噪比為?25 dB，結合圖11 的誤比特率曲線，能夠到達接近 10?3的誤比特率性能。MUSIC 算法與CBF 算法性能接近（假設MUSIC 算法的信源數(shù)已知），這2 種算法最低工作信噪比分別為?17 dB和?18 dB，波束成形的導向正確時可以取得10 dB的增益，也就是說將MUSIC 算法或者CBF 算法應用到擴頻通信中，信噪比為?17 dB 和?18 dB 左右，可以實現(xiàn)接近 10?3的誤比特率性能。ISM 和CSM-TCT 抗噪聲性能較差，最低工作信噪比分別為?13 dB 和?8 dB，對比圖11，將ISM 或者CSM-TCT 應用于擴頻通信中可能會導致擴頻通信性能的損失。

3.4 抗多普勒分析

由于水聲通信會受到發(fā)射機和接收機的相對運動的影響，導致信號壓縮或者擴展，故利用重采樣模擬多普勒效應，設置壓縮因子a=0.002 5。同時，假設水聲擴頻信號的入射角為30°，此時將信噪比設置較高，觀察多普勒帶來的影響，SNR=0，其他條件保持不變。

多普勒效應對空時分簇DoA 估計影響較大，若圖2 所示的算法框架中沒有加時頻二維搜索，直接將多波束成形后的數(shù)據(jù)做相關，得到的空時分簇DoA 估計如圖17 所示，與圖12 沒有多普勒效應的空時分簇DoA 估計對比可知，多普勒效應不僅會導致峰值的幅度下降，而且會使估計結果變得混亂。

圖17 無時頻二維搜索的空時分簇DoA 估計

根據(jù)時頻二維搜索的原理，7 階擴頻信號的頻率分辨率約為9.8 Hz，壓縮因子a=0.002 5時，頻偏約為10 Hz，此時擴頻信號相關輸出會變小且沒有明顯峰值。若按照圖2 所示的算法框架加入時頻二維搜索空時分簇DoA 估計會出現(xiàn)明顯峰值，且空時信道結構更加清晰，如圖18 所示。

圖18 有時頻二維搜索的空時分簇DoA 估計

綜上所述，本文所提空時分簇DoA 估計算法利用m 序列的偽隨機特性，對干擾信號具有天然抑制作用，并將方向譜上升至時延?角度譜，使分辨率更加高，利用m 序列優(yōu)良的相關性，估計均方根誤差較低，具有較高的估計精度和穩(wěn)定性，空時分簇DoA 估計方法繼承了常規(guī)波束成形的優(yōu)良特性，結構簡單，估計相干信號源不受影響。由于直接對接收水聲擴頻信號進行空時二維搜索會受到多普勒效應的影響，故加入時頻二維搜索，解決多普勒效應的影響后可以穩(wěn)定工作。總之，相比與CBF 算法、文獻[11]中的MUSIC 算法、ISM、CSM-TCT，基于水聲擴頻信號的空時分簇DoA 估計方法性能更加優(yōu)越。

4 實驗驗證

為了驗證本文算法的有效性，在哈爾濱工程大學水聲實驗室進行水池實驗，發(fā)射擴頻信號中心頻率為4 kHz，帶寬為2.5 kHz，擴頻序列為9 階m序列，采用一個發(fā)射換能器發(fā)射信號，陣元間距為18 cm 的10 元均勻直線陣列接收信號，實驗設備如圖19 所示。均勻直線陣列垂直放置在水深5 m 的信道水池中，陣列中心位于水下2 m 左右，發(fā)射換能器與陣列中心等高，兩者相距20 m。

圖19 實驗設備

實測水池的空時分簇DoA 估計如圖20 所示。信號以1°入射到均勻直線陣上，經過水面水底的反射，?17°與26°有一些信號簇到達，圖20 也體現(xiàn)了水池的空間信道結構。

圖20 實測水池的空時分簇DoA 估計

實測水池的信道估計如圖21 所示，以1°波束成形后的信道估計如圖22 所示。對比圖21 和圖22 可知，波束成形后信號幅度增加10 倍，且削弱了其他方向的多徑干擾。同時，空時分簇DoA 估計得到的信道結構是三維的，可以估計出某條路徑的入射方向。因此，空時二維搜索可以得到更清晰的信道結構，為通信系統(tǒng)開發(fā)創(chuàng)造更多的先驗知識。

圖21 實測水池的信道估計

圖22 以1°波束成形后的信道估計

根據(jù)圖23 可知，由于信道水池中多徑干擾嚴重，對相干信號源敏感的ISM 已經失效。CBF 算法的空間分辨率由于“瑞利限”的限制，空間譜主瓣寬度明顯增加，但是仍然能觀察到信號主要從?17°、1°和26°入射。CSM-TCT 初始預估角度采用CBF 估計得到的角度，分辨率有明顯的改善。由于信源數(shù)未知，CSM-TCT 和MUSIC 算法可能會出現(xiàn)虛假峰值。本文所提空時分簇DoA 估計算法更穩(wěn)健，且基本上每個入射角均可以在圖20(a)所示的空時分簇三維圖上找到。

圖23 水池信號的對比算法DoA 估計結果

5 結束語

本文提出了一種基于水聲擴頻信號的空時分簇DoA 估計算法。該算法思想取自時頻二維搜索，實質上是將空間所有方向進行波束搜索，利用擴頻序列的自相關性，與各個方向波束成形后的數(shù)據(jù)做相關，將方向譜上升至時延?角度譜，使其保持較高的分辨率。由于所提算法利用了m 序列的偽隨機特性，因此不受干擾信號的影響。擴頻帶來的擴頻增益使估計均方根誤差較低，在低信噪比下具有較高的估計精度和穩(wěn)健性。本文首先對原理進行推導分析，然后將空時分簇DoA 估計算法與CBF 算法、文獻[11]中的MUSIC 算法、ISM、CSM-TCT 進行對比，仿真結果表明，空時分簇DoA 估計算法具有更加穩(wěn)健的性能?？諘r分簇DoA 估計算法可提升水聲擴頻信號DoA 估計中的抗干擾性、分辨率、抗噪聲以及穩(wěn)健的性能，既保留了擴頻的各種優(yōu)點，又兼具了陣列體制的優(yōu)勢，但其存在計算量較大的問題，期望未來得以解決。