潘浩, 孫超, 卓頡, 劉雄厚

(西北工業(yè)大學(xué) 聲學(xué)工程研究所, 陜西 西安 710072)

水下聲成像系統(tǒng)(如扇掃聲吶系統(tǒng)、側(cè)掃聲吶系統(tǒng)和合成孔徑聲吶系統(tǒng)等)已廣泛應(yīng)用于水下目標(biāo)探測(cè)、海底地形測(cè)繪和海洋礦物勘探等一系列軍用和民用領(lǐng)域。為了提高成像系統(tǒng)的方位向分辨率,一些高分辨成像聲吶系統(tǒng)通常采用由成千上萬(wàn)個(gè)陣元組成的平面陣或柱面陣,從而導(dǎo)致其成本和運(yùn)算量居高不下[1]。于是,對(duì)成像陣列進(jìn)行稀疏優(yōu)化研究便成為聲吶成像領(lǐng)域的熱點(diǎn)之一。文獻(xiàn)[2-4]采用模擬退火(simulated annealing,SA)算法、文獻(xiàn)[5-6]采用遺傳算法(genetic algorithm,GA)研究了傳統(tǒng)陣列(線(xiàn)列陣和平面陣)的稀疏優(yōu)化方法。然而,為了滿(mǎn)足陣列波束圖主瓣寬度和旁瓣級(jí)的要求,所稀疏的陣元數(shù)目往往是很有限的,對(duì)于陣元基數(shù)本身就特別大的傳統(tǒng)高分辨成像聲吶系統(tǒng)來(lái)說(shuō),稀疏化后的系統(tǒng)成本和運(yùn)算量依然不夠理想[7]。

為了獲得更好的節(jié)省陣元效果,MIMO技術(shù)被引入到陣列成像領(lǐng)域[7-10]。文獻(xiàn)[8-9]建立一種由相互垂直的2條均勻線(xiàn)列陣組成的MIMO陣列模型進(jìn)行雷達(dá)成像,得出了與同等成像分辨力的矩形平面接收陣相比,該MIMO陣列可以節(jié)省M×N-(M+N)個(gè)陣元(其中M為發(fā)射陣元個(gè)數(shù),N為接收陣元個(gè)數(shù))的結(jié)論。鑒于聲吶系統(tǒng)與雷達(dá)系統(tǒng)在信號(hào)處理方法上有許多相似之處,已有學(xué)者將MIMO陣列應(yīng)用于水下聲成像領(lǐng)域,文獻(xiàn)[7]在陣列尺寸上優(yōu)化了文獻(xiàn)[8-9]的MIMO陣列模型,建立了一種由4條均勻線(xiàn)列陣圍成的口字型MIMO聲吶模型,并用來(lái)進(jìn)行水下三維成像;并通過(guò)設(shè)計(jì)多種用于水下三維成像的MIMO陣列布陣和推導(dǎo)MIMO虛擬陣元的解析坐標(biāo)。文獻(xiàn)[10]也將MIMO技術(shù)引入聲吶成像領(lǐng)域,利用近場(chǎng)聚焦寬帶波束形成成像算法對(duì)MIMO成像聲吶進(jìn)行了仿真和水池試驗(yàn),驗(yàn)證了MIMO聲吶可以提高方位分辨率和抗混響能力。

與傳統(tǒng)的成像聲吶系統(tǒng)相比,MIMO成像聲吶系統(tǒng)利用波形分集技術(shù)產(chǎn)生大量的虛擬陣元,從而用較少的實(shí)際陣元個(gè)數(shù)和較小的物理陣列孔徑就能夠達(dá)到與前者同樣的成像效果[7]。對(duì)于造價(jià)成本昂貴并且易受安裝平臺(tái)限制的水下聲成像系統(tǒng)來(lái)說(shuō),MIMO聲吶有很大的發(fā)展空間。然而,MIMO成像聲吶系統(tǒng)在獲得虛擬陣元優(yōu)勢(shì)的同時(shí)意味著接收端匹配濾波運(yùn)算量的增加,影響了系統(tǒng)的實(shí)時(shí)成像性能?；诖?將稀疏陣列技術(shù)應(yīng)用到MIMO聲吶二維成像當(dāng)中,提出了利用稀疏陣列技術(shù)來(lái)降低MIMO聲吶二維成像中匹配濾波處理運(yùn)算量的方法。針對(duì)MIMO聲吶虛擬接收陣列為矩形平面陣的情況,采用2條相互垂直的均勻線(xiàn)列陣(uniform linear array,ULA)構(gòu)成的MIMO陣,使用模擬退火算法對(duì)該虛擬接收陣的陣元位置和加權(quán)系數(shù)同時(shí)優(yōu)化,使獲得的波束圖逼近預(yù)先給定的期望波束圖,同時(shí)最小化開(kāi)啟的虛擬陣元數(shù)目,進(jìn)而移除關(guān)閉的虛擬陣元對(duì)應(yīng)的匹配濾波器。此時(shí)MIMO聲吶實(shí)際物理陣列保持不變,但匹配濾波處理的通道數(shù)大幅度減少。這樣就從系統(tǒng)根源上改善了MIMO二維成像聲吶系統(tǒng)的性能,加快了成像聲吶的數(shù)據(jù)處理速率和實(shí)時(shí)成像速率。

1 MIMO聲吶布陣方式及其虛擬接收陣列稀疏化

1.1 MIMO聲吶布陣方式

本文采用文獻(xiàn)[8-9]提到的2條均勻線(xiàn)列陣互相垂直組成的MIMO聲吶。假設(shè)有M個(gè)發(fā)射陣元位于x軸上,N個(gè)接收陣元位于y軸上,發(fā)射陣列和接收陣列的陣元間距均為λ/2(其中λ為發(fā)射信號(hào)的中心頻率在水下聲速1 500 m/s時(shí)對(duì)應(yīng)的波長(zhǎng))且都以坐標(biāo)原點(diǎn)為中點(diǎn)。MIMO聲吶的三維坐標(biāo)系統(tǒng)如圖1所示,其中實(shí)心圓代表發(fā)射陣元,空心圓代表接收陣元,θ表示俯仰角,φ表示方位角。

由文獻(xiàn)[7]可知:遠(yuǎn)場(chǎng)條件下,MIMO陣列的虛擬發(fā)射陣元位于原點(diǎn),虛擬接收陣元坐標(biāo)等于一對(duì)實(shí)際的發(fā)射、接收陣元坐標(biāo)之和。因此,本文采用的MIMO聲吶布陣的虛擬接收陣列為M×N的均勻矩形平面陣,如圖2所示,其中圖2a)為MIMO聲吶實(shí)際物理陣列,圖2b)為其虛擬陣列。

圖1 MIMO聲吶布陣示意圖 圖2 MIMO聲吶物理陣列及其虛擬陣列布陣形式

1.2 MIMO聲吶虛擬接收陣稀疏化

陣列的稀疏優(yōu)化是指在一個(gè)均勻滿(mǎn)陣中關(guān)掉一些陣元和改變一些陣元的加權(quán)系數(shù),在波束圖滿(mǎn)足成像所需分辨率和加權(quán)系數(shù)比(current taper ratio,CTR)(CTR定義為陣元最大加權(quán)系數(shù)與最小加權(quán)系數(shù)的比值)閾值條件的前提下,求得此時(shí)開(kāi)啟陣元數(shù)目的最小值[3]。其中,加權(quán)系數(shù)比的設(shè)定能減輕具有最大幅度加權(quán)系數(shù)的陣元發(fā)生問(wèn)題而導(dǎo)致的影響,保證穩(wěn)健性。這里定義坐標(biāo)系為u-v空間:

u=sinθcosφ

(1)

v=sinθsinφ

(2)

式中:θ和φ分別為俯仰角和方位角。

采用陣元位置和幅度加權(quán)同時(shí)優(yōu)化的方法,定義模擬退火算法的能量函數(shù)如(3)式所示[4]

(3)

式中:X是陣元的位置矩陣,W是幅度加權(quán)系數(shù)矩陣,p(u,v)/Q表示當(dāng)前歸一化的波束圖,bd為一個(gè)常數(shù),表示期望的波束圖的最大旁瓣級(jí)。S為在期望旁瓣區(qū)域內(nèi)滿(mǎn)足不等式p(u,v)/Q>bd所有(u,v)值的集合,Q為所有幅度加權(quán)系數(shù)wi的總和,A為當(dāng)前開(kāi)啟陣元的個(gè)數(shù),Ro表示稀疏優(yōu)化過(guò)程中得到的當(dāng)前加權(quán)系數(shù)比(obtainedCTR),Rd表示期望目標(biāo)加權(quán)系數(shù)比(desiredCTR),k1、k2和k3是常數(shù),決定各參量的權(quán)重大小。優(yōu)化的目的就是通過(guò)迭代找到合適的幅度加權(quán)矩陣W,使得當(dāng)前的波束圖p(u,v)的旁瓣區(qū)域逼近期望的最大旁瓣級(jí)bd、Ro逼近Rd的同時(shí)得到開(kāi)啟的最少陣元數(shù)目Amin,最終的優(yōu)化結(jié)果W是這3個(gè)參數(shù)的折中。

當(dāng)波束主瓣掃描時(shí),每個(gè)位置上的陣元按相位差公式(4)補(bǔ)償相位[3]

(4)

式中：ri表示第i個(gè)陣元的位置坐標(biāo),a0表示波束圖主瓣指向的方向向量。因此,第i個(gè)陣元對(duì)應(yīng)的復(fù)加權(quán)系數(shù)wi可表示為

wi=wi·exp(-j·φi),i=1,2…,MN

(5)

式中：wi為幅度加權(quán),φi為相位加權(quán)。

圖3所示為采用模擬退火算法對(duì)矩形平面陣進(jìn)行稀疏優(yōu)化的流程圖。

圖3 模擬退火算法流程圖

2 低運(yùn)算量MIMO聲吶二維成像方法

設(shè)遠(yuǎn)場(chǎng)中的待成像區(qū)域可建模為P個(gè)離散點(diǎn),MIMO聲吶發(fā)射正交多相編碼信號(hào),那么第n個(gè)接收陣元上的接收信號(hào)可以表示為[7]

(6)

在接收端,每個(gè)接收陣元都使用M個(gè)匹配濾波器分別對(duì)M個(gè)發(fā)射波形進(jìn)行匹配,從而分離出這M個(gè)發(fā)射信號(hào)。這樣一來(lái),就需要使用MN個(gè)匹配濾波器,而一個(gè)匹配濾波輸出相當(dāng)于產(chǎn)生一個(gè)虛擬陣元。因此,MIMO成像聲吶系統(tǒng)在獲得虛擬陣元優(yōu)勢(shì)的同時(shí)帶來(lái)了接收端匹配濾波處理的大運(yùn)算量。圖4所示為匹配濾波器與虛擬陣元一一對(duì)應(yīng)關(guān)系,其中空心圓代表實(shí)際接收陣元,陰影圓代表虛擬接收陣元。根據(jù)這樣的一一對(duì)應(yīng)特點(diǎn),結(jié)合上一節(jié)虛擬接收陣列的稀疏化結(jié)果,移除關(guān)閉虛擬陣元對(duì)應(yīng)的匹配濾波器。匹配濾波器數(shù)目的減少意味著匹配濾波處理運(yùn)算量和參加波束形成算法的原始數(shù)據(jù)的減少,MIMO成像聲吶系統(tǒng)的運(yùn)算量因而就可以從本質(zhì)上得到降低。

圖4 虛擬接收陣元與匹配濾波器一一對(duì)應(yīng)關(guān)系

(7)

式中：n0代表保留下來(lái)的虛擬陣元編號(hào),xtm和xrn分別為該虛擬陣元對(duì)應(yīng)的第m個(gè)實(shí)際發(fā)射陣元和第n個(gè)實(shí)際接收陣元坐標(biāo)。由于匹配濾波器與虛擬接收陣元是一一對(duì)應(yīng)的關(guān)系,根據(jù)(8)式可知,MIMO聲吶中的第n個(gè)接收陣元上的第m個(gè)匹配濾波器需要開(kāi)啟。

設(shè)虛擬陣列稀疏后的陣元個(gè)數(shù)為N0。根據(jù)(8)式,第n0(n0=1,2,…,N0)個(gè)匹配濾波輸出yn0(t),可表示為

yn0(t)=xn(t)*hm(t)

(8)

式中：*代表卷積,hm(t)為與第m個(gè)發(fā)射信號(hào)對(duì)應(yīng)的匹配濾波器的沖擊響應(yīng)函數(shù),其表達(dá)式為

hm′(t)=[sm(T-t)]H

(9)

式中：[]H為取共軛,T=LT0為單個(gè)發(fā)射信號(hào)的長(zhǎng)度(T0為發(fā)射信號(hào)單個(gè)子碼長(zhǎng)度,L為子碼個(gè)數(shù))。

由以上分析可知,經(jīng)過(guò)稀疏優(yōu)化后,MIMO聲吶的虛擬接收陣元個(gè)數(shù)從MN減少為N0,需要的匹配濾波器數(shù)量也從MN減少為N0。對(duì)(9)式中的N0個(gè)匹配濾波輸出進(jìn)行移相波束形成,那么第q(q=1,2,…,Q)個(gè)波束輸出Bq(t)的表達(dá)式為

(10)

式中：wn0為第n0個(gè)虛擬陣元對(duì)應(yīng)的復(fù)加權(quán)值。

最后調(diào)節(jié)陣列指向使得目標(biāo)區(qū)域被多個(gè)窄波束覆蓋,按照回波到達(dá)的先后順序提取各個(gè)波束下的強(qiáng)度,獲得目標(biāo)區(qū)域的二維強(qiáng)度圖。

3 仿真實(shí)驗(yàn)

3.1 虛擬接收陣的稀疏優(yōu)化仿真

采用24發(fā)、24收的MIMO陣列和圖1給出的布陣方式,并利用模擬退火算法對(duì)其虛擬陣列(24×24的均勻平面陣)進(jìn)行優(yōu)化。根據(jù)圖1所示的MIMO聲吶布陣方式和探測(cè)區(qū)域,優(yōu)化后的三維波束圖俯仰角θ的范圍為[90°,180°],方位角φ的范圍為[0°,360°],坐標(biāo)變換到u-v空間,預(yù)設(shè)期望波束圖的最高旁瓣級(jí)bd為0.1(20*log10(0.1)=-20 dB),期望旁瓣范圍滿(mǎn)足u2+v2>0.08,期望權(quán)重系數(shù)比Rd=6.0,其他相關(guān)設(shè)置如下:Tstart=1 000,k1=1 000,k2=0.2,k3=2 000,同時(shí)優(yōu)化陣元位置和幅度加權(quán)系數(shù),將576個(gè)陣元稀疏到263個(gè)。優(yōu)化前后的陣型如圖5所示,為了說(shuō)明稀疏后的陣元位置在空間多個(gè)波束指向上均可以保證其性能,圖6分別給出了稀疏優(yōu)化后平面陣的主瓣(θ0,φ0)分別指向(120°,60°)、(150°,120°)和(180°,150°)時(shí)多主瓣指向波束圖,對(duì)比起見(jiàn),一并給出稀疏前均勻陣列的側(cè)視圖。圖7給出優(yōu)化后保留陣元對(duì)應(yīng)的加權(quán)系數(shù)幅度值的三維圖。

圖5 平面陣列稀疏前后的陣元位置

圖6 多主瓣指向三維波束圖

圖7 幅度加權(quán)系數(shù)的三維圖

3.2 低運(yùn)算量MIMO聲吶二維成像仿真

MIMO陣列的發(fā)射信號(hào)為多相編碼信號(hào),其子碼個(gè)數(shù)L=64,單個(gè)子碼長(zhǎng)度T0=100Ts(Ts=1/fs且fs=500 kHz),載頻f0=100 kHz。假設(shè)水下有8個(gè)理想散射點(diǎn)位于z=-5 m的平面上且呈“F”型,位置坐標(biāo)分別為(-1,1,-5)m、(-1,0,-5)m、(-1,-1,-5)m、(-1,-2,-5)m、(0, 1,-5)m、(1, 1,-5)m、(1,0,-5)m和(0,0,-5)m,如圖8所示。

圖8 散射點(diǎn)與成像陣列的相對(duì)位置

其中圖8a)中“+”號(hào)代表MIMO陣列所處的位置。散射點(diǎn)散射強(qiáng)度均設(shè)為1。在接收端,接收陣元上的信噪比設(shè)為10 dB,噪聲為加性高斯白噪聲。波束俯仰角從120°掃到180°,方位角從0°掃到360°,間隔均為3°,共形成2 541個(gè)波束。二維成像結(jié)果如圖9a)所示。為了對(duì)比起見(jiàn),給出未稀疏虛擬陣列(即24×24的均勻平面陣)的MIMO聲吶二維成像的仿真結(jié)果如圖9b)所示,其中發(fā)射信號(hào)與散射點(diǎn)位置同上,幅度加權(quán)采用-20 dB等旁瓣級(jí)加權(quán)。

圖9 MIMO聲吶二維成像結(jié)果

從圖9a)、圖9b)對(duì)比看出,采用本文提出的低運(yùn)算量MIMO聲吶成像方法和未稀疏虛擬陣的常規(guī)MIMO聲吶成像方法都能夠清楚地顯示散射點(diǎn)的位置,但由于利用了模擬退火算法將虛擬陣元個(gè)數(shù)從576稀疏到263,相應(yīng)地所須匹配濾波器的個(gè)數(shù)由576減少為263個(gè),運(yùn)算量降低了大約一半,比較2個(gè)方法的運(yùn)行時(shí)間也可以看出來(lái)這一點(diǎn)。仿真工具為Matlab2011b版軟件,程序運(yùn)行的硬件平臺(tái)為雙核3.30 GHz主頻、4G內(nèi)存PC機(jī),為保證計(jì)時(shí)準(zhǔn)確,只對(duì)運(yùn)算時(shí)間較長(zhǎng)的匹配濾波處理部分進(jìn)行計(jì)時(shí),并作了20次重復(fù)試驗(yàn),耗時(shí)結(jié)果見(jiàn)圖10。

圖10 重復(fù)試驗(yàn)20次的耗時(shí)曲線(xiàn)

可見(jiàn),對(duì)于24發(fā)、24收的MIMO聲吶來(lái)說(shuō),采用未稀疏虛擬陣MIMO聲吶成像方法耗時(shí)大約33 s,而低運(yùn)算量MIMO聲吶成像方法僅耗時(shí)大約15 s。

4 結(jié) 論

高成本與大運(yùn)算量一直是成像聲吶系統(tǒng)面臨的2大難題。本文提出的MIMO聲吶低運(yùn)算量二維成像方法,既利用了MIMO陣列節(jié)省陣元的能力以降低成像系統(tǒng)的造價(jià)成本,又運(yùn)用稀疏陣列技術(shù)在保證成像分辨能力的前提下有效減少了匹配濾波器的數(shù)目。該方法從根源上降低MIMO聲吶成像系統(tǒng)的運(yùn)算量,提高了成像聲吶的數(shù)據(jù)處理速率和實(shí)時(shí)成像速率。由于稀疏陣列技術(shù)發(fā)展已較為成熟,該方法理論上可以應(yīng)用到任意布陣方式的MIMO陣列二維和三維成像系統(tǒng)中,如文獻(xiàn)[10]提到的由1條發(fā)射直線(xiàn)陣與1條接收弧形陣組成的MIMO陣列(其虛擬陣列為柱面陣),以及由2條弧形陣(一條弧形陣為發(fā)射陣,另一條弧形陣為接收陣)組成的MIMO陣列(其虛擬陣列為類(lèi)似球面陣的曲面陣)等,具有一定的通用性。

參考文獻(xiàn)：

[1] Murino V, Trucco A. Three-Dimensional Image Generation and Processing in Underwater Acoustic Vision [J]. Proc IEEE, 2000, 88(12): 1903-1948

[2] Trucco A, Murino V. Stochastic Optimization of Linear Sparse Arrays [J]. IEEE Journal of Oceanic Engineering, 1999, 24(3): 291-299

[3] Trucco A. Thinning and Weighting of Large Planar Arrays by Simulated Annealing [J]. IEEE Trans on Ultrasonic, Ferroelectrics and Frequency Control, 1999, 46(2):347-355

[4] Chen P, Zheng Y Y, Zhu W. Optimized Simulated Annealing Algorithm for Thinning and Weighting Large Planar Arrays in both Far-Field and Near-Field [J]. IEEE Journal of Oceanic Engineering, 2011, 36(4): 658-664

[5] O’Neill D J. Element Placement in Thinned Arrays Using Genetic Algorithms [C]∥Proceedings of IEEE OCEANS’94. Brest: IEEE, 1994: 301-306

[6] Haupt R L. Thinned Arrays Using Genetic Algorithms [J]. IEEE Trans on Antennas and Propagation, 1994, 42(7): 993-999

[7] Liu X H, Sun C, Yi F, et al. Underwater Three-Dimensional Imaging Using Narrowband MIMO Array [J]. Science China in Physics Mechanics and Astronomy, 2013, 56(7): 1346-1354

[8] Wang H J, Su Y. Narrowband MIMO Radar Imaging with Two Orthogonal Linear T/R Arrays [C]∥International Conference on Signal Processing, Beijing, 2008: 2513-2516

[9] Wang D W, Ma X Y, Su Y. Two-Dimensional Imaging via a Narrowband MIMO Radar System with Two Perpendicular Linear Arrays [J]. IEEE Trans on Image Processing, 2010, 19(5): 1269-1279

[10] 滕婷婷,孫大軍,劉鑫,田原. 波形分集MIMO成像聲吶技術(shù)研究[J]. 哈爾濱工程大學(xué)學(xué)報(bào), 2013, 34(5): 581-587

Teng T T, Sun D J, Liu X, Tian Y. The Techniques for Waveform Diversity MIMO Imaging Sonar[J]. Journal of Marine Science and Application, 2013, 34(5): 581-587 (in Chinese)