基于主動時反的淺海目標DOA估計優化算法

荊海霞, 王海燕, 劉鄭國, 申曉紅, 張之琛

1.西北工業大學 航海學院, 陜西 西安 710072; 2.西安外事學院 工學院, 陜西 西安 710077;3.中國船舶重工集團公司 辦公廳, 北京 100097

波達方向(direction of arrival, DOA)估計是陣列信號處理的一個重要研究課題，在雷達、聲吶、通信等領域有著廣泛的應用。一般來說，DOA估計方法大都是建立在直達波模型上，忽略多徑或想辦法消除多徑帶來的影響[1-3]，其把多徑看作是一個不利的因素；在直達波模型基礎上，提出了各種各樣的波達方向估計方法，但這些方法大都建立在高信噪比條件下，在低信噪比時算法性能下降，甚至失效。

近些年提出的時間反轉[4](time reversal，TR或時反)方法，由于利用了聲互易性及時反不變性原理，可自適應地修正多徑信道引起的信號畸變，使目標在源位置獲得空時聚焦，因此可將時反方法應用于目標的DOA估計中，嘗試解決低信噪比情況下常規方法無法正確估計的問題。

其中，文獻[5-6]將主動時反(active time reversal,ATR)方法分別應用到探地雷達和移動通信網絡的目標DOA估計中，各自推導了TR/DOA估計器的CRB，并和未引入TR的DOA估計器性能進行了比較；仿真結果表明，引入TR方法后DOA估計器的性能得到了提高。文獻[7]在文獻[6]的基礎上，加入了距離估計器，實現了對單目標的主動定位。文獻[8]在文獻[7]的基礎上，研究了多輸入多輸出(MIMO)系統中的目標DOA問題，建立了TR/MIMO模型；隨后文獻[9]又將TR引入多普勒頻移條件下的MIMO雷達目標DOA估計中；文獻[10]將壓縮感知(compressive sensing,CS)和TR技術結合應用到MIMO雷達，實現了DOA、DOD(direction of departure)及多普勒頻移的聯合估計。文獻[11]提出了一種虛擬時間反轉(virtual time reversal,VTR)DOA估計算法，應用于空間電磁輻射源的遠場窄帶信號源方位角的被動估計中，通過查找掃描區域中能量最大的點來確定輻射源的方位角。以上電磁波領域的時反DOA估計方法，對實現水下目標的時反DOA估計有一定的借鑒意義。

文獻[12]針對均勻淺海目標DOA估計，提出了一種基于非均勻線列陣(non-uniform line array,NLA)的被動時反(passive time reversal,PTR)超指向性模型，從信號檢測角度建立了仿真所需模型，采用常規波束形成方法，實現了低信噪比情況下目標的方位估計。文獻[13]針對水下聲多徑時變信道，提出了一種多普勒效應下的TR目標定向算法，通過頻率補償與時間反轉處理后，可以精確地實現目標定向。以上水下目標時反定向方法建立在被動時反方法基礎上，首先需解決如何精確獲取真實海洋信道的問題，因為只有當模擬信道與實際海洋信道完全匹配時，才能達到理想的TR空時聚焦效果，從而達到精確DOA估計的目的。

主動時間反轉相對被動時間反轉，由于利用了海洋信道傳播的互易性，無需環境的先驗知識，可簡化信號處理的運算量，相比而言更適用于惡劣環境，物理上更容易實現。

基于此，本文將結合前期研究的主動時反聚焦特性[14-15]和探測特性[16]，進一步利用主動時間反轉方法對水下目標的俯仰角DOA估計問題進行研究，基于水聲信道的射線理論和陣列信號處理中的直達徑模型，建立了基于均勻線列陣(uniform line array,ULA)的常規多徑DOA估計和主動時反DOA估計模型，應用Capon算法對所建模型進行了DOA估計實現，對有無主動時反時的估計算法性能進行了比較。仿真結果表明，在同樣的信噪比情況下，所提出的主動時反DOA方法可以獲得比常規方法更好的估計結果，尤其是在低信噪比情況下更為明顯。

1 均勻線列陣多徑DOA估計模型

基于射線理論的均勻線列陣多徑DOA估計如圖1所示。

圖1 均勻線列陣多徑DOA估計模型

圖中，收發合置換能器陣(source-receive array,SRA)為均勻線列陣，其陣元個數為P，陣元間距為d；為方便表示，圖中只畫出了目標與SRA之間的3條傳播路徑：直達波、海面反射波與海底發射波，其入射角分為θ,α,β。但為更具有一般性，后面推導的聲傳播路徑不止局限于這3條徑，而采用N條徑表示。

1.1 常規多徑DOA估計模型

位于k陣元的PS(探測聲源)發射信號f(t),根據射線理論,設發射聲源與目標之間的信道傳輸函數為:

hk(t)=∑Nn=1cknδ(t-τkn)

(1)

式中,N表示聲線總數,ckn,τkn分別表示第k個陣元和目標之間第n條本征聲線(也可認為是傳播路徑)對應的衰減幅度和時延。忽略探測聲源到目標的信道噪聲,則目標接收到的信號為:

xk(t)=f(t)?hk(t)=∑Nn=1cknf(t-τkn)

(2)

忽略目標反射的影響(即假定目標反射系數為1),考慮接收過程噪聲的影響,則SRA第j個陣元接收的信號為:

yj(t)=xk(t)?hj(t)+vj(t)=

∑Nn=1cknf(t-τkn)?∑Mm=1cjmδ(t-τjm)+vj(t)=

∑Mm=1∑Nn=1cjmcknf(t-τjm-τkn)+vj(t)

(3)

式中,hj(t)表示目標與第j個陣元之間的信道函數,vj(t)表示第j個陣元接收的環境噪聲,cjm,τjm分別表示目標和第j個陣元之間的第m條路徑對應的衰減幅度和時延。根據對海洋聲場的仿真可知,目標經第m條徑到達各個陣元的衰減幅度差別非常小,為了簡化表示,可認為其只和路徑有關,而與陣元無關,即可將cjm表示成cm。

設發射信號形式為:f(t)=s(t)ejωct,結合陣列信號處理的遠場窄帶模型理論,(3)式可表示為:

yj(t)=∑Mm=1∑Nn=1cmcnf(t-τkn-τ1m-Δτjm)+vj(t)≈

∑Mm=1∑Nn=1cmcns(t)ejωc(t-τkn-τ1m-Δτjm)+vj(t)=

∑Mm=1∑Nn=1cmcns(t)ejωcte-jωcτkne-jωcτ1me-jωcΔτjm+vj(t)=

∑Mm=1cme-jωcΔτjme-jωcτ1m·∑Nn=1cne-jωcτkn·s(t)ejωct+vj(t)

(4)

式中,τ1m表示的是目標到達陣元1的第m條徑的時延,Δτjm表示的是目標經過第m條徑到達j陣元與第1陣元的相對時延。

將第j個陣元接收信號推廣到各個陣元,并表示成矩陣形式,(4)式可變為:

Y(t)=ACDXkF(t)+V(t)

(5)

式中,Y(t)=[y1(t),…,yP(t)]T表示SRA的1-P個陣元接收到的信號;A為P×M矩陣,表示相對首陣元的時延矩陣,如(6)式所示;C為M階對角矩陣,M為目標到各陣元的路徑總數,其對角元素cm表示目標經過第m條路徑到達各陣元的衰減幅度;F(t)=s(t)ejωct稱為發射矩陣;D=[e-jωcτ11,e-jωcτ12,…,e-jωcτ1M]T表示的是目標到陣列1的時延矩陣。

Xk=∑Nn=1cne-jωcτkn=A(k,:)CD

可認為是從k陣元發射的信號經過信道到達目標的接收信號矩陣;V(t)表示噪聲矩陣。

A=e-jωcΔτ11e-jωcΔτ12…e-jωcΔτ1M

e-jωcΔτ21e-jωcΔτ22…e-jωcΔτ2M

????

e-jωcΔτP1e-jωcΔτP2…e-jωcΔτPM

(6)

參考陣列信號處理理論,當均勻線列陣的陣元間距為d時,可將(6)式表示成:

A=[a(θ1),a(θ2),…,a(θM)]=

11…1

e-j2πλdsinθ1e-j2πλdsinθ2…e-j2πλdsinθM

????

e-j2πλ(P-1)dsinθ1e-j2πλ(P-1)dsinθ2…e-j2πλ(P-1)dsinθM

(7)

A為包含了多徑信息的陣列流形矩陣(又稱:方向矩陣),主要取決于陣列結構與目標來波方向,其中第m列各量代表的是目標經過第n條徑到達各陣元的信息;θM表示目標到陣元的第m條徑的角度信息。

1.2 基于ATR的多徑DOA估計模型

考慮主動時反情況,將SRA各陣元接收信號分別進行時反,以第j個陣元為例,有:

yj(-t)=B·s(-t)e-jωct+vj(-t)

(8)

式中,B=∑Mm=1cme-jωcΔτjme-jωcτ1m·∑Nn=1cne-jωcτkn,其值與時間t無關,可看成一系數;又由于采用主動時反,發射信號頻率已知,故可采用濾波器處理[16],在時反前先消除噪聲的影響,故(8)式可變為:

yj(-t)≈B·s(-t)e-jωct=Bf(-t)

(9)

將其作為二次發射信號重新發射到信道中,依然滿足遠場窄帶模型理論;重復上述第一次過程,參考(1)～(4)式,設時反后第l個陣元接收信號為zl(t),有:

zl(t)≈∑Mm=1∑Nn=1cmcnBs(-t)e-jωc(t-τjn-τ1m-Δτlm)+

vl(t)=∑Mm=1cmejωcΔτlmejωcτ1m·∑Nn=1cnejωcτjn·

Bs(-t)e-jωct+vl(t)

(10)

再參考(5)～(7)式,可得第j個陣元時反發射后SRA接收信號Zj(t)為:

(11)

式中,符號“*”表示共軛,C,A,D,V同(5)式,Xj和(5)式中的Xk意義相同,只不過Xk為k陣元發射的接收信號矩陣,Xj為j陣元發射的接收信號矩陣。

將SRA的第j個陣元時反發射情況推廣到所有陣元,可得時反后SRA的接收信號總和為:

(12)

2 DOA估計算法

2.1 常規Capon算法

根據陣列信號處理理論,Capon算法的目的是想辦法減小噪聲和期望信號以外的干擾信號功率,同時把期望方向的信號功率維持在一定增益上[17],結合(5)式建立的模型可將Capon算法表示成以下最小值問題:

minwH(θ)YCw(θ)

s.t.:wH(θ)a(θ1)=1

(13)

利用Lagrange乘子法可求得上述問題的解為:

將(14)式帶入(13)式可得，期望方向以外功率最小時期望方向的功率最大，為：

P=1aH(θ1)θ1)(15)

因此Capon算法的空間譜可定義為：

P=1aH(θ)(θ)(16)

通過對(16)式進行譜峰搜索可求出目標的DOA值。式中：

a(θ)=[1,e-j2π/λdsinθ,…,e-j2π/λ(P-1)dsinθ]T(17)

2.2 主動時反Capon算法

根據(12)式，遵循上述常規DOA估計同樣的步驟，有：

(θ)Ztrwtr(θ)

Ptr=1aH(θ)(θ)(19)

通過對(19)式進行譜峰搜索可求出目標的DOA值。

相對常規Capon算法，主動時反算法中的目標反射信號，由于在接收陣進行了時反操作后重新發射到目標上，根據時反的聚焦特性可知，SRA各陣元的二次發射信號將會在目標處形成聚焦，該過程相當于波束形成過程將波束聚焦在了目標上，只不過時反方法由于利用了多徑，相對波束形成其聚焦在目標上的能量更大。同時，文獻[16]論證了單陣元主動時反探測陣可提高接收信號的信噪比，此結論可推廣至多陣元探測。從以上2個角度可知，利用主動時反方法進行目標的DOA估計，其估計的精度會更高，抑制旁瓣的能力會更強。

3 仿真研究

考慮均勻淺海波導環境，各深度聲速均為1.5 km/s。仿真模型參考圖1，其中時反陣(SRA)取12個陣元，陣元間距為0.75 m，首陣元1#距離水面75 m；PS的深度為78.75 m，目標深度165.5 m，PS與目標的水平距離為1 km，海深300 m；PS發射頻率為1 kHz 的CW信號，快拍數取1 000。

按以上條件，在MATLAB仿真環境下，利用Bellhop專用仿真工具箱模擬海洋聲場環境，忽略海面或海底的多次反射路徑，只考慮圖1所示的3條傳播路徑情況，可得到目標與SRA各陣元之間的各路徑的時延、幅值及角度信息。通過仿真可知，目標反射的聲波經海面一次反射徑(或直達徑，或海底一次反射徑)到達各陣元的幅值和角度信息非常接近，可近似將同一路徑的值取為一個值；按以上路徑順序，到達各陣元的幅度信息依次約為：{0.65×10-3,1.0×10-3,0.5×10-3}，角度依次約為：{13.68°, -5.00°, -19.58°}；目標到參考陣元(首陣元)的時延值依次為：{0.686, 0.669, 0.708}，其他參數略。

為了對比同樣條件下引入主動時反前后Capon算法的性能，信噪比定義為：考慮多徑時的首次接收信號中信號與噪聲的比值；主動時反二次接收信號加入和首次接收信號一樣的噪聲。

圖2和圖3分別為信噪比分別為-15 dB和-20 dB時常規Capon算法和加入主動時反后的ATR Capon算法的DOA估計圖，圖中豎線表示角度期望值。

圖2 SNR=-15 dB時Capon算法和ATR Capon 算法的DOA估計圖

圖3 SNR=-20 dB時Capon算法和ATR Capon 算法的DOA估計圖

從圖2可以看出，在信噪比為-15dB時，ATR Capon算法的旁瓣遠低于其相應主瓣能量，估計值更接近于待估計目標期望值，且分辨率高于常規Capon算法；當信噪比為-20 dB時，常規Capon算法已不能正確估計出目標，而ATR Capon算法估計出的目標值基本不變。

為了更好地比較出以上2種算法的性能，對上述2種算法分別做了1 000次蒙特卡羅仿真，比較在不同信噪比情況下，2種算法估計的均方根誤差(RMSE)情況，如圖4所示。從圖4中可以看出，當信噪比大約從-12 dB開始2種算法的RMSE非常接近，故圖5放大了-12～0 dB的情況。

圖4 Capon算法和ATR Capon算法的 DOA估計均方根誤差分析圖

綜合圖4和圖5可以看出，ATR Capon算法估計出的角度均方根誤差要小于常規Capon算法，尤其是當信噪比非常低時，ATR Capon算法的優勢尤為明顯。

圖5 信噪比為-12～0 dB時Capon算法和ATR Capon 算法的DOA估計均方根誤差分析圖

4 結 論

本文利用主動時間反轉方法研究了淺海目標DOA估計性能。仿真結果表明：在低信噪比情況下，基于ATR Capon方法的DOA估計相對于常規Capon方法的DOA，分辨率更高，抑制旁瓣的能力更強，估計結果更準確，可很好地應用于低信噪比多徑條件下的DOA估計。

參考文獻：

[1] Daeipour E, Blair W D, Bar-Shalom Y. Bias Compensation and Tracking with Monopulse Radars in the Presence of Multipath[J]. IEEE Trans on Aerospace and Electronic Systems, 1997, 33(3): 863-882

[2] Garber S M. High Resolution Sonar Signals in a Multipath Environment[J]. IEEE Trans Aerospace and Electronic Systems, 1966, AES-2(6): 431-440

[3] Derryberry J H, Gregg W D. On Optimizing Array Reception of Multipath[J]. IEEE Trans on Aerospace Electronic System, 1970, AES-6(2): 188-199

[4] Mathias Fink. Time Reversal of Ultrasonic Fields-Part I: Basic Principles[J]. IEEE Trans on Ultrasonics Ferroelectrics and Frequency Control, 1992, 39(5): 555-566

[5] Foroohar Foroozan, Amir Asif. Cramer-Rao Lower Bound for Time Reversal Active Array Direction of Arrival Estimation in Multipath Environment[C]∥2010 IEEE International Conference on Acoustics, Speech and Signal Prolessing: 2646-2649

[6] Foroohar Foroozan, Amir Asif. Time Reversal Direction of Arrival Estimation with Cramer-Rao Bound Analysis[C]∥2010 IEEE Global Telecommunications Conference: 1-5

[7] Foroohar Foroozan, Amir Asif. Time Reversal Based Active Array Source localization[J]. IEEE Trans on Signal Processing, 2011, 59(6): 2655-2668

[8] Foroohar Foroozan, Amir Asif. Direction Finding Algorithms for Time Reversal MIMO Radars[C]∥2011 IEEE Statisticol Signal Processing Workshop: 433-436

[9] Foroohar Foroozan, Amir Asif. Time Reversal MIMO Radar for Angle-Doppler Estimation[C]∥2012 IEEE Statistical Signal Processing Workshop: 860-863

[10] Mohammad H S, Sajjadieh, Amir Asif. Compressive Sensing Time Reversal MIMO Radar: Direction and Doppler Frequency Estimation[J]. IEEE Tran on Signal Processing Letters, 2015, 22(9): 1283-1287

[11] Fu Yongqing, Liu Wei, Bai Ruijie, et al. A Novel Virtual Time Reversal Method for Passive Direction of Arrival Estimation[J]. Mathematical Problems in Engineering, 2015(2015):1-12

[12] 楊伏洲,王海燕,申曉紅,等. 基于時間反轉的非均勻線列陣超指向性陣元分布模型[J]. 上海交通大學學報, 2013,47(12): 1907-1910

Yang Fuzhou, Wang Haiyan, Shen Xiaohong, et al. Super-Direction Element Distribution Model of NLA Based on TR[J]. Journal of Shanghai Jiaotong University, 2013, 47(12): 1907-1910 (in Chinese)

[13] Shao Jianfeng, Zhang Xiaomin, Liu Yihai, et al. Estimation of Time Reversal Target DOA over Underwater Acoustic Multipath Time-Varying Channel[C]∥2014 IEEE China Summit and International Conference on Signal and Infomation Processing: 795-799

[14] 荊海霞,李洪義. 基于主動時間反轉的水下目標自適應聚焦研究[J]. 電子設計工程,2015,23(24):12-15

Jing Haixia, Li Hongyi. Self-Adaptive Focusing of Underwater Target Based on Active Time Reversal[J]. Electronic Design Engineering, 2015, 23(24):12-15 (in Chinese)

[15] 荊海霞,申曉紅,王海燕. 淺海聲源信道中信號優化探測研究[J]. 計算機仿真,2016,33(10):157-161

Jing Haixia, Shen Xiaohong, Wang Haiyan. Signal Optimization Detection in Acoustic Source Channel under Shallow Water[J]. Computer Simulation, 2016, 33(10): 157-161 (in Chinese)

[16] 荊海霞,申曉紅,劉鐳. 基于主動時間反轉的目標探測性能研究[J]. 電視技術, 2016,40(8): 103-107

Jing Haixia, Shen Xiaohong, Liu Lei. Performance Study of Target Detection Based on Active Time Reversal[J]. Video Engineering, 2016,40(8): 103-107 (in Chinese)

[17] 孫超. 水下多傳感器陣列信號處理[M]. 西安：西北工業大學出版社，2007

Sun Chao. Underwater Multi-sensor Array Signal Processing[M]. Xi′an, Northwestern Polytechnical University Press, 2007 (in Chinese)