張羽， 王朋， 劉紀元， 鐘榮興， 韋琳哲， 遲騁

(1.中國科學院 聲學研究所， 北京 100190； 2.中國科學院 先進水下信息技術重點實驗室， 北京 100190；3.中國科學院大學， 北京 100049)

0 引言

合成孔徑聲納(SAS)是一種高分辨率成像聲納，其原理是利用小孔徑基陣沿方位向勻速直線運動、形成一個虛擬的大孔徑，對虛擬大孔徑內的回波解耦進行匹配濾波處理，從而得到不隨距離變化而僅與基陣尺寸相關的方位向高分辨率。SAS實現(xiàn)高分辨率成像的關鍵是平臺沿著理想直線航跡做勻速運動，但在實際工作中，海浪、暗涌或自身運動不穩(wěn)定等因素會影響聲納平臺的航行軌跡，導致其產生運動誤差。運動誤差會造成聲納成像結果出現(xiàn)畸變、散焦、重影等現(xiàn)象，降低圖像質量。因此，需要對運動誤差進行補償[1]。運動補償按照運動誤差數(shù)據(jù)源的不同，分為基于運動測量系統(tǒng)的運動補償、基于回波數(shù)據(jù)的運動補償和圖像域自聚焦?；谶\動測量系統(tǒng)的運動補償是利用安裝在聲納平臺的高精度運動傳感器測量并記錄聲納的姿態(tài)、速度等信息，通過坐標轉換、積分等運算獲得聲納的運動軌跡，從而得到運動誤差。其中，運動測量系統(tǒng)通常包括慣性導航系統(tǒng)(INS)、多普勒測速儀(DVL)、全球定位系統(tǒng)(GPS)和深度計等?；诨夭〝?shù)據(jù)的運動補償是對前后兩幀回波信號做相關運算，利用相關系數(shù)反演出運動誤差，稱為基于重疊相位中心(DPC)算法。圖像域自聚焦則是在圖像域提取相位誤差并將其消除[2-3]。

近年來，國內的SAS運動補償研究主要沿著運動測量系統(tǒng)和回波數(shù)據(jù)兩個方面分別展開。文獻[4-5]分別研究了改進型粒子濾波和卡爾曼濾波在SAS運動測量系統(tǒng)中的應用，通過多傳感器數(shù)據(jù)融合方法提高運動誤差的估計精度，該方法難點在于選取合適的系統(tǒng)模型。文獻[6]考慮了回波沿斜距向的空變性，對傳統(tǒng)基于回波數(shù)據(jù)的DPC方法進行改進，采用分段DPC擬合出運動誤差，對大測繪帶成像聲納的改善效果較好。文獻[7]針對雙側成像聲納，提出了基于雙側回波DPC方法，區(qū)分出不同類型的運動誤差，獲得了比單側補償更好的成像結果。

國外的SAS運動補償研究則以運動測量系統(tǒng)與DPC方法相結合為主。文獻[8-9]將DPC看作輔助傳感器，對其數(shù)據(jù)預處理后輸入濾波模型，與INS、DVL、GPS、深度計等多種運動傳感器數(shù)據(jù)進行融合估計運動誤差，甚至在某些水聲環(huán)境良好的時刻，直接用DPC取代運動傳感器計算運動誤差。文獻[10-11]認為DPC結果中包含著海底地形起伏和介質傳播特性等信息，是實現(xiàn)精確補償必須要考慮的關鍵要素。因此，使用DPC數(shù)據(jù)與INS數(shù)據(jù)之間的差值修正INS，得到了更準確的運動誤差估計結果。文獻[12]基于拖曳式低頻SAS開展了DPC、INS與GPS 3種傳感器數(shù)據(jù)融合的運動補償研究，其中GPS數(shù)據(jù)來自水面艦。雖然在外參標定過程中采用了較多的近似處理，融合過程略顯粗糙，但是成像結果依然表明了數(shù)據(jù)融合的有效性。

此外，國內外一些機構還開展了DPC聯(lián)合INS服務于導航系統(tǒng)的研究，其目的是將DPC技術應用在水下導航領域。雖然應用背景有別于運動補償，但二者技術手段非常相似，都采用多傳感器數(shù)據(jù)融合的方法，利用DPC估計橫蕩誤差和縱蕩誤差精度高的特點，彌補INS長時間工作容易產生積累誤差的缺點。文獻[13]研究了搭載SAS的水下無人平臺在無GPS情況下的自主導航問題，使用DPC輔助INS獲得了比DVL輔助INS更準確的航行軌跡。類似地，文獻[14]比較了DVL、DPC、陀螺儀、INS等多種傳感器不同組合方式對導航結果造成的影響，研究結果驗證了DPC輔助INS的可行性和優(yōu)越性。

運動補償?shù)暮脡呐c運動誤差估計的準確性直接相關，無論是運動測量系統(tǒng)還是DPC，本質上都是為了獲得聲納的實際航跡，從而計算實際航跡與理想航跡的偏差。隨著對聲納回波研究程度的不斷加深，越來越多的DPC相關算法被相繼提出，魯棒性和準確性也在逐漸提高[15]。文獻[16-17]充分利用多子陣優(yōu)勢，使用回波數(shù)據(jù)空間互相關矩陣估計聲納前行速度，可以達到與DVL測速相近的結果。另外，考慮到運動測量系統(tǒng)往往包含多種傳感器，系統(tǒng)復雜度高且存在數(shù)據(jù)冗余。如果聲納自身數(shù)據(jù)能夠輔助測量運動參數(shù)，代替功能重復的傳感器，則不僅有利于降低系統(tǒng)復雜度，還可以節(jié)省空間、節(jié)約成本。

基于以上背景，本文提出一種SAS DPC+INS運動誤差估計算法，首先使用DPC聲相關測速取代DVL多普勒頻移測速的核心功能，獲得聲納的前行速度，然后根據(jù)矢量運算法則換算出聲納在載體坐標系下的三向速度，最后建立卡爾曼濾波模型，融合DPC的速度與INS的加速度，獲得聲納速度的最優(yōu)估計值，用于計算運動誤差。通過分析試驗數(shù)據(jù)發(fā)現(xiàn)，DPC+INS方法可以有效提高運動誤差估計準確度，改善目標重影、散焦的情況，提高圖像質量。

1 多子陣SAS聲相關測速

為了提高測繪效率，SAS通常采用多子陣的配置，即單個發(fā)射陣和多個接收子陣搭配的模式。多子陣的配置促使SAS領域產生了利用回波數(shù)據(jù)的相關性估計聲納運動參數(shù)的DPC方法。沿距離向對等效相位中心的兩幀回波數(shù)據(jù)分析，重疊相位中心的相關函數(shù)峰值到相關序列中心的偏移量反映了前后兩幀信號之間的時延，換算成空間的距離就是聲納的橫蕩誤差。沿方位向對等效相位中心的兩幀回波數(shù)據(jù)分析，通過相關函數(shù)最大值確定重疊的等效相位中心的位置，根據(jù)等效相位中心之間的間隔計算出聲納沿方位向移動的距離，即可得到縱蕩誤差。本文使用的DPC方法屬于后者，對回波信號沿方位向處理獲取縱蕩誤差，進而換算成聲納前向速度。

根據(jù)波形不變性理論[18]，只要SAS相鄰兩幀的等效相位中心存在重疊，通過對這兩幀回波數(shù)據(jù)求相關，利用相關系數(shù)即可解析出聲納的運動參數(shù)[19-20]。在多子陣空間互相關矩陣中，重疊相位中心對的個數(shù)與脈沖重復周期以及SAS基陣的前向速度，三者之間存在等式關系。通過等式關系得到聲納基陣的前向速度后，按照速度矢量分解法則，結合偏航角和俯仰角，便可計算出聲納前向(方位向)、距離向和垂向3個方向上的速度。

1.1 等效相位中心假設

SAS一般采用單發(fā)多收的配置，其中發(fā)射陣與接收陣是分置的。為了處理方便，引入“等效相位中心”假設，把分置的發(fā)射陣和接收陣等效成一個收發(fā)共置的虛擬陣元，該虛擬陣元的位置處于發(fā)射陣和接收陣的中心，即假定聲波的發(fā)射路徑和接收路徑是一致的。如圖1(a)所示。圖1中，A點為發(fā)射陣T和接收陣R的等效相位中心，rt為真實的發(fā)射路徑，rr為真實的接收路徑，聲程r′的2倍近似為rt與rr之和，d為發(fā)射陣和接收陣之間的間隔，B為目標。由等效相位中心假設所引入的聲程誤差與d和r′有關。對于收發(fā)間距分米量級，探測距離百米量級的成像聲納，聲程誤差一般較小，可忽略不計[21]。圖1(b)給出了單發(fā)多收陣的一系列等效相位中心示意圖。

圖1 SAS等效相位中心Fig.1 SAS equivalent phase center

1.2 多子陣空間互相關矩陣估計前向速度

選取合適的前行速度，保證SAS等效相位中心在前后兩幀之間存在重疊，即可利用回波數(shù)據(jù)的相關性解算出聲納斜距向和方位向的位移。如圖2所示，假設聲納存在N個等距排列的等效相位中心，編號依次為{1，2，3，…，N-1,N}，相鄰兩個相位中心之間的間距為l. 聲納以vy(y表示方位向)的速度沿方位向勻速前進，從第k-1幀時刻到第k幀時刻，聲納沿方位向產生了位移Δy，由于運動誤差的存在，聲納同時在斜距向上產生了位移Δx′(x′表示斜距向)。

圖2 重疊相位中心聲相關測速原理Fig.2 Principle of acoustic correlation velocity measurement by DPC

假定在這兩幀時刻之間有p對相位中心完全重疊，定義PRT為第k-1幀與第k幀的時間間隔，稱為脈沖重復周期，則重疊相位中心對的個數(shù)、脈沖重復周期以及聲納前進速度存在如下等式關系：

Δy=(N-p)·l=vy·PRT，

(1)

因此，只要確定了等效相位中心重疊的個數(shù)，就可以計算出聲納的前向速度。

根據(jù)波形不變性理論，如果前后兩幀有重疊的等效相位中心，則這兩個相位中心的回波數(shù)據(jù)具有較大的相關性。按照(2)式對第k-1幀第i個通道的回波信號與第k幀第j個通道的回波信號取互相關，得到互相關系數(shù)ci,j，

(2)

式中：1≤i≤N；1≤j≤N；corr[·]表示相關運算；si,k-1(m)表示第k-1幀通道i接收到的回波信號，m為距離向時間變量；sj,k(n)表示第k幀通道j接收到的回波信號，n為距離向時間變量。

(3)

圖3 脈壓后回波信號互相關矩陣Fig.3 Spatial mutual correlation matrix of echoes after pulse-compression

必須指出，利用回波數(shù)據(jù)空間互相關矩陣僅能估計出聲納的前向速度。當聲納在前后兩幀出現(xiàn)如圖2所示的不平行情況時，該方法等效為把第k幀時刻基陣的位置(實線)近似為與前一時刻平行的位置(虛線)。

1.3 聲納三向速度關系式

實際工作中，聲納在前后幀之間不僅會產生方位向的位移，而且受水流或自身運動不穩(wěn)定等因素的影響，聲納還會產生斜距向的位移，該位移由橫蕩誤差和升沉誤差組成，分別沿著聲納載體坐標系的距離向和垂向。將一對重疊的等效相位中心得到的兩組數(shù)據(jù)進行相關，通過相關峰與相關序列中心之間的偏移量可計算出斜距向位移，但傳統(tǒng)的DPC方法并不能夠直接區(qū)分出橫蕩誤差和升沉誤差。文獻[6]采用分段DPC和最小二乘法擬合出橫蕩誤差和升沉誤差，文獻[7]采用雙側DPC方法通過解方程的形式計算出橫蕩誤差和升沉誤差，均增加了計算量?？紤]到DPC方法本身存在一定的估計誤差，多次運算疊加容易產生積累誤差，而SAS配備有高精度INS，能夠輸出角度誤差數(shù)據(jù)，結合1.2節(jié)計算出的前向速度，對速度矢量按照矢量運算法則進行分解，利用前向速度和角度誤差解算出距離向速度和垂向速度。

聲納與INS之間屬于剛性連接，運動過程中具有相同的角度誤差，包括偏航角、俯仰角、橫滾角。下面將聲納看成1個點，在載體坐標系下對其速度矢量v進行分析。如圖4所示，聲納前進方向指向y軸，φh、φp分別表示偏航角與俯仰角，可將速度矢量v分解為沿距離向x、方位向y和垂向z3個方向的速度，分別為vx、vy、vz，其與v之間的等式關系可表示為

(4)

圖4 聲納速度矢量分解圖Fig.4 Vector decomposition diagram of sonar velocity

由此，根據(jù)1.2節(jié)計算的前向速度vy可以得到

(5)

2 DPC+INS算法

準確地，vx、vy、vz是DPC方法估計出來的速度，相比于聲納真實速度仍存在誤差，而INS可以輸出聲納的加速度。按照多傳感器數(shù)據(jù)融合的思想，采用卡爾曼濾波技術對二者的數(shù)據(jù)進行融合，能夠克服單個傳感器的不確定性，提高速度數(shù)據(jù)估計的準確度，稱為DPC+INS算法。對速度積分得到聲納的實際航跡，采用最小二乘法擬合出理想直線航跡，運動誤差即為實際航跡與理想航跡之間的偏差。

2.1 DPC+INS模型

圖5 DPC與INS數(shù)據(jù)融合模型Fig.5 Data fusion model of DPC+INS

DPC聯(lián)合INS的目的是將DPC估計的方位向、距離向、垂向速度與INS測量的方位向、距離向、垂向加速度進行融合，圖5給出了DPC+INS數(shù)據(jù)融合模型。首先對回波數(shù)據(jù)進行相關運算，估計出聲納的方位向速度vy；然后結合INS輸出的姿態(tài)角對速度矢量分解，解算出聲納的距離向速度vx和垂向速度vz；接著將距離向、方位向、垂向三向速度連同載體坐標系下INS輸出的距離向、方位向、垂向加速度ax、ay、az以及姿態(tài)角一并輸入卡爾曼濾波器進行濾波；最后輸出聲納速度估計值，同樣包括距離向、方位向和垂向3個方向的數(shù)據(jù)。

2.2 DPC+INS數(shù)據(jù)融合方法

在數(shù)據(jù)融合前需要對數(shù)據(jù)進行預處理，由于載體坐標系是隨著平臺的運動而變化的，為了方便后續(xù)計算運動誤差，需要使用固定不變的坐標系。因此，按照(6)式將載體坐標系下的速度數(shù)據(jù)vx、vy、vz和加速度數(shù)據(jù)ax、ay、az分別與姿態(tài)轉換矩陣相乘，得到地理坐標系下的速度數(shù)據(jù)vn、ve、vu和加速度數(shù)據(jù)an、ae、au：

(6)

式中：vn、ve、vu分別表示北向速度、東向速度和天向速度；an、ae、au分別表示北向加速度、東向加速度和天向加速度；Th、Tp、Tr分別為偏航角φh、俯仰角φp和橫滾角φr的轉換矩陣，角度誤差數(shù)據(jù)來自INS.

根據(jù)速度與加速度之間的關系，DPC+INS系統(tǒng)方程可描述為線性方程，量測量為DPC的三向速度和INS的三向加速度，濾波時間間隔為1個脈沖重復周期PRT.

使用完全狀態(tài)法建立狀態(tài)方程如下：

(7)

式中：

(8)

(9)

(10)

建立量測方程如下：

(11)

式中：Z6×1表示vn、ve、vu、an、ae、au的量測值；H6×6表示6×6階的量測系數(shù)矩陣，

(12)

η6×1表示量測噪聲，為白噪聲。

圖6 卡爾曼濾波流程Fig.6 Kalman filtering process

根據(jù)線性系統(tǒng)理論，把狀態(tài)方程和量測方程離散化，選定部分參數(shù)初始值?？柭鼮V波遵循最小均方誤差準則對量測值進行遞推濾波，按照圖6所示的流程迭代，不斷更新狀態(tài)估計值。一個濾波周期包括了時間更新和量測更新兩個過程，分別體現(xiàn)了狀態(tài)預測和狀態(tài)估計的思想。圖6中：0為狀態(tài)量初始估計值；P0為估計均方誤差初始值；下標k為采樣時刻；Pk為估計均方誤差；Φk|k-1為一步轉移陣，由F6×6離散化得到；Γk-1為系統(tǒng)噪聲驅動陣；Qk和k分別為系統(tǒng)噪聲方差陣及估計值；Pk|k-1為一步預測均方誤差；k為狀態(tài)量6×1離散化的結果；k|k-1為狀態(tài)一步預測值；Zk為量測量Z6×1離散化后的結果；Hk為量測陣；Rk和k分別為量測噪聲方差陣及估計值；Kk為濾波增益；I為單位矩陣。

經過卡爾曼濾波后得到狀態(tài)量的估計值，即為DPC數(shù)據(jù)與INS數(shù)據(jù)融合后的結果。選取狀態(tài)量的前3個變量{n,e,u}作為聲納東北天三向的速度，相當于把DPC作為主數(shù)據(jù)源，而INS的作用是在DPC的基礎上不斷增補和豐富一些信息。

2.3 運動誤差的計算與補償

多子陣SAS的運動誤差分為角度誤差和平動誤差，對其中單個陣元而言，角度誤差本質上是以平動誤差的形式影響回波數(shù)據(jù)。因此，運動補償?shù)年P鍵是計算出平動誤差。根據(jù)文獻[6]的分析，平動誤差中橫蕩誤差和升沉誤差對回波相位影響最大，相位誤差Δψ近似表示為

(13)

式中：λ為聲波波長；Δx為橫蕩誤差；Δh為升沉誤差；α為聲納掠射角。

可見相位誤差是關于橫蕩誤差和升沉誤差的函數(shù)，下面將對實際航跡分析，計算橫蕩誤差和升沉誤差。如圖7所示，設初始時刻的位置O為原點，對聲納東北天三向速度{n,e,u}積分可得到聲納的實際航跡{sn,se,su}，記為曲線s，第k幀時刻的實際位置為Sk；利用最小二乘擬合方法，獲得聲納理想的直線航跡{s′n,s′e,s′u}，記為s′，第k幀時刻的理想位置為S′k；目標在B點. 則根據(jù)圖 7虛線框內的幾何關系，橫蕩誤差和升沉誤差分別為

Δx=|SkB|cosα-|S′kB|cosα，

(14)

Δh=|SkB|sinα-|S′kB|sinα.

(15)

圖7 聲納運動軌跡與誤差模型Fig.7 Model of sonar trajectory and motion error

將橫蕩誤差和升沉誤差代入(13)式，得到相位誤差，通過頻域補相的方式補償原始回波，即可實現(xiàn)運動誤差的補償。

3 數(shù)據(jù)分析

為驗證DPC+INS算法的有效性，選取2018年項目組在某水庫典型地貌區(qū)域采集的高頻SAS數(shù)據(jù)進行分析。聲納系統(tǒng)參數(shù)如表1所示，根據(jù)表1中參數(shù)，結合第2節(jié)的知識可知，等效相位中心間距為0.05 m. 實際工作時，前后兩幀的相位中心約有一半可重疊，為前向速度的估計提供了重要支撐。多子陣空間互相關矩陣結果參考圖3，不難發(fā)現(xiàn)，前一幀通道1～通道11與后一幀通道14～通道24形成1條亮度最大的直線，表明該條直線上的相關系數(shù)最大。

表1 SAS系統(tǒng)參數(shù)Tab.1 System parameters of SAS

圖8給出了該航次的姿態(tài)角變化曲線。由圖8可以看出，整段數(shù)據(jù)在80 s時間內航向角變化約15°，屬于變化較大的情況，而俯仰角和橫滾角的變化較小且相對平穩(wěn)。

圖8 某航次聲納姿態(tài)角變化曲線Fig.8 Changing curves of sonar attitude angle

圖9 聲納三向速度曲線Fig.9 Curves of sonar three-dimensional velocity

圖9給出了使用DPC+INS算法得到的北向、東向和天向速度估計結果。為了便于對比，引入DPC方法形成對照。從圖9中的曲線趨勢可以看出，DPC+INS算法估計的結果與DPC方法估計的結果大體趨勢相同。細致觀察圖9(a)和圖9(b)不難發(fā)現(xiàn)，DPC速度曲線比較平滑，而DPC+INS相比于DPC，曲線上帶有許多毛刺，這是因為DPC的速度融合了INS加速度，增加了很多細節(jié)成分。圖9(c)的天向速度曲線中，DPC+INS算法振蕩明顯較大，經過分析，發(fā)現(xiàn)是INS天向通道數(shù)據(jù)不穩(wěn)定引起的。但由于天向速度數(shù)值較小，天向位移(升沉誤差)對回波的影響并不如東北向位移(橫蕩誤差)大。

圖10 地貌成像結果Fig.10 Topography imaging results

使用多子陣ωK算法[22-23]對該航次數(shù)據(jù)進行成像，結果如圖10所示。其中圖10(a)是未經運動補償?shù)脑汲上駡D，受運動誤差的影響。從圖10(a)可以看到：框選區(qū)域1和2處目標均存在方位向的重影，單個目標成像結果卻是2個，相當于引入了虛假目標；框選區(qū)域3的目標同時存在散焦和重影的現(xiàn)象。圖10(b)為使用DPC方法補償?shù)膱D像，相比于圖10(a)，解決了框選區(qū)域1和2方位向重影的問題，但框選區(qū)域3的底部依然存在較嚴重的散焦。圖10(c)是使用DPC+INS算法補償?shù)膱D像，整體成像效果較好，與圖10(a)相比，解決了框選區(qū)域1和2處的目標方位向重影的問題；與圖10(b)相比，框選區(qū)域3處目標散焦的情況也有所改善，表明DPC結合INS算法估計運動誤差更準確。也就是說，將INS數(shù)據(jù)與DPC數(shù)據(jù)融合，可以達到比單獨使用DPC算法更好的效果。

圖11 目標2方位向剖面圖Fig.11 Azimuth profile of Target 2

為了客觀地評價不同算法的補償效果，下面對目標的成像質量做進一步比較。參照合成孔徑成像性能指標，給出圖10中目標2的方位向剖面圖，結果如圖11所示。從圖11(a)中可以看出，未經運動補償?shù)哪繕?在方位向上存在多個波瓣，能量較分散，且副峰的峰值較高，因此出現(xiàn)了圖11(a)中目標重影的情況。觀察圖11(b)和圖11(c)可以發(fā)現(xiàn)，經過運動補償后，目標2的能量集中在主瓣，聚焦效果得到明顯改善。相比之下，DPC+INS算法的補償效果比DPC算法更好一些，不僅降低了峰值旁瓣比，還減小了主瓣寬度。

表2給出了不同情況下目標2方位向的性能指標，包括峰值旁瓣比和主瓣寬度。通過對比性能指標可知，DPC+INS算法運動補償效果的最好。

表2 目標2方位向性能指標Tab.2 Azimuth performance index of Target 2

4 結論

為了降低系統(tǒng)復雜度，提高SAS回波數(shù)據(jù)的利用率，本文將DPC方法與INS相結合，提出一種DPC+INS算法。利用多子陣空間互相關矩陣解算出聲納前向速度，通過INS姿態(tài)角計算出DPC三向速度，并采用卡爾曼濾波融合DPC三向速度與INS三向加速度，求解聲納速度的最優(yōu)值計算運動誤差?；诤嚁?shù)據(jù)對算法進行了驗證。所得主要結論如下：

1) DPC結合INS估計三向速度是可行的，DPC+INS算法相比于DPC方法，融合INS的加速度數(shù)據(jù)，增加細節(jié)信息，提高了聲納速度估計的準確性，為運動誤差的計算提供了重要支撐。

2) 通過對比原始數(shù)據(jù)成像圖，DPC+INS算法和DPC方法均能提高圖像質量，解決目標重影的問題。相比之下，DPC+INS算法補償效果更好，能夠有效改善目標散焦的情況。