劉曉，李海森，周天，徐超，么彬

(哈爾濱工程大學水聲技術重點實驗室，黑龍江哈爾濱150001)

多波束測深聲吶作為一種高效、高精度的海洋勘測設備，具有可同時測深和海底成像的優(yōu)點.特別是近年來其成像結果在海底生境調查、底質分類等研究中都得到了廣泛的應用［1-2］.目前，利用多波束測深聲吶對海底成像主要有3種方法:1)形成2個額外的寬波束，對回波信號進行類似側掃聲吶的高密度幅度時間序列采樣，稱為偽側掃成像(pseudosidescan imagery)［3］;2)每個接收窄波束只取一個聲強值，如波束主軸方向的聲強值［4］，或者是每個波束內的平均聲強值［5］;3)對每個接收窄波束都進行幅度時間序列采樣，稱為“腳印時間序列(footprint time series)”或“片斷(snippet)”法，其與前2種方法相比，避免了聲強數據與測深數據的融合問題，能同時獲得高信噪比與高分辨率的聲吶圖像，因此應用更為廣泛［6］.然而，由于傳統(tǒng)的多波束測深算法，如幅度或相位檢測法只能估計出每個波束主軸方向海底檢測點處的地理位置和時延信息，所以在“Snippet”方法中每個波束內除主軸方向其他聲強樣本的地理位置是通過假設波束內為平海底情況得到的，而這種理想假設使得這些聲強值與地理信息不能一一對應，特別是在波束照射區(qū)域內的地形起伏較大時，如海溝、海山等，成像質量可能受到較大影響.由于多子陣檢測法在海底地形探測時能獲得較高精度以及良好的空間分辨能力［7］，所以進一步考慮將其應用到海底成像技術中.利用該方法獲得大量海底檢測點的回波到達時間(TOA)和到達方向(DOA)信息，并由此估計各檢測點的空間位置與回波強度，目的是使生成的海底圖像在具有較高分辨能力的前提下，實現海底地貌與地形的共點測量.

1 海底成像方法

1.1 多子陣檢測法原理

多子陣檢測法是一種相位檢測法，主要思想是將一個多元均勻直線接收陣形成多個相互有一定重疊的子陣，對每個子陣分別進行多波束形成，計算各子陣的相位函數，由此估計出相位斜率，并通過圖像變換算法得到海底各檢測點處回波的TOA和DOA.由于文獻［7］對此方法有較為詳細的論述，所以這里只對其基本原理作簡要說明.

各個子陣經波束形成后，能夠提供每個波束內回波信號的相位信息，稱為子陣相位函數，一般表示為

式中:φk(θ(n)，θr)為第 k個子陣的相位函數;K(n，θr)為相位斜率，簡記為K，它是時間n和波束控制角θr的函數;xk為第k個子陣等效聲中心的坐標位置，β為第1個陣元接收信號的相位.對所有子陣的相位信息進行最小二乘擬合運算就能估計出K.

根據式(1)可求出相鄰2個子陣之間的相位差，結合由基陣理論所得的相位差函數，便可得到信號到達角序列θ(n)，若忽略基陣安裝傾角，則θ(n)可表示為

利用海底圖像變換算法將相位斜率圖像K(n，θr)變換為時間和信號到達角的函數B(n，θ)，圖像中的像素點值表征了海底某方向回波的概率.變換前，將B(n，θ)中所有的像素值都初始化為零，然后利用式(2)對K(n，θr)中的每一像素值計算到達角θ，從而可獲得每個采樣時刻對應的回波DOA估計值，有效提高了其空間分辨率.此外，圖像變換時對鄰近波束的統(tǒng)計平均處理能有效提高DOA的估計精度.

1.2 空間位置解算

在聲速已知情況下，利用1.1節(jié)中得到的TOA和DOA估計值可進一步解算出各檢測點的海底深度和水平位移.但由于聲波穿透水體到達海底經歷的各個水層中的聲速間存在差異，這樣導致了實際的聲波傳播路線為曲線而非直線，因此海底檢測點處的空間位置并不能用三角關系進行簡單的計算，而需要在聲速剖面基礎上，采用層追加的方法來解算各檢測點相對于船體坐標系的坐標位置［8］.

如圖1所示，假設波束經歷的水柱有N+1層，聲速剖面內有 N+1 個聲速采樣點(C0，1，…，N)，層 i的上下表面聲速分別為 Ci－1和 Ci(i=1，2，…，N)，該層厚度為Δzi.為保證計算精度，這里采用基于層內常梯度假設下的聲線跟蹤技術，即假設聲速在層i內以常梯度gi變化，可表示為

式中:gi為常數.則第i層內聲線的水平位移為

式中:p由Snell法則得到:

式中:θ0為估計所得的DOA值，C0為聲基陣表面聲速.

圖1 層內常梯度聲線跟蹤示意Fig.1 Schematic diagram of sound velocity correction based on layered constant velocity gradient model

聲線在該層經歷的時間ti為

其中CHi為Harmonic平均聲速.

在實際的計算中，將傳播時間ti逐層累加并與實際測量的單程傳播時間，即TOA值的一半，進行比較判斷傳播層數，從而解算出各海底檢測點相對于基陣中心的水平位移;再結合基陣的GPS數據(如經度、緯度)與測量船的姿態(tài)信息(如艏向、橫搖、縱搖、升沉)，就可以估計出這些檢測點的空間位置坐標，即地形信息.由于海底各檢測點是利用多子陣檢測法得到的實際測點，所以避免了“Snippet”方法中當波束內的地形存在起伏時解算的檢測點位置與實際情況存在誤差的問題.

1.3 回波強度估計

對于多波束測深聲吶而言，將其接收基陣做多波束形成，并在每個接收窄波束輸出信號的基礎上估計海底各檢測點處的回波強度可以提高成像數據的信噪比以及空間分辨能力.考慮到基陣波束指向性對接收信號幅度的影響，只保留－3 dB波束寬度范圍內的波束輸出信號用于回波強度的估計［9］，如圖2所示.回波持續(xù)時間與水深、波束控制角以及波束寬度有關.對于θr方向的波束，－3 dB功率點對應的波束上限角度θh和下限角度θl分別為

式中:N是陣元長度，d是陣元間距，λ是波長.由此可得－3 dB波束寬度為

假設在θr方向上，接收波束在其形成的海底波束腳印內共有M個檢測點被多子陣檢測法檢測到，且其 TOA 和 DOA 估 計 值 分 別 為 T1，2，…，i，…，M和θ1，2，…，i，…，M，如圖 2.利用這些 TOA 值對 － 3 dB 波束寬度內的反向散射信號進行采樣，得到M個長度等于發(fā)射脈寬大小的幅度時間序列，并由此估計各檢測點對應的回波強度.如，第r號波束內第i個檢測點處的回波強度可表示為

式中:E(n)是第i個時間序列的幅度，τ是發(fā)射脈沖寬度，fs是采樣頻率，n是采樣樣本數，［ ］表示取整.

圖2 成像數據測量原理示意Fig.2 Schematic diagram of imaging datamea surement

根據上文可知，該方法進行回波強度估計所利用的DOA和TOA與進行空間位置解算時所用的是同一序列，所以這兩者的測量針對的是相同的海底檢測點，實現了共點測量.

1.4 成像數據的修正

聲波在傳播過程中受到傳播損失TL、波束指向性D、有效的聲照射區(qū)域面積A以及系統(tǒng)處理增益G等因素影響，根據式(5)計算得到的回波強度I并不能直接反映海底地貌特征，因此在實際的成像數據測量中需對它們進行補償以獲得海底反向散射強度BS，根據聲吶方程可得

式中:I0是參考聲強，SL是聲源級.

引起聲強在介質中傳播衰減的主要原因包括擴展損失和吸收損失，由于聲波波陣面按球面擴展，所以總傳播損失等于:

式中:r為傳播距離，α為聲吸收系數.

由于每個波束只保留－3 dB波束寬度內的樣本強度，所以只對式(3)和式(4)角度范圍內的海底檢測點進行波束圖修正.由乘積定理可知歸一化的基陣波束指向性為

每個波束接收到的信號來自波束腳印內的聲照射區(qū)域，對于內側波束，波束腳印大小小于聲照射區(qū)域，但隨著入射角度的增加，波束腳印會不斷增大并逐漸超過聲照射區(qū)域，如圖3所示，所以A應取聲照射區(qū)域面積Ainsonif和波束腳印面積Abeam的相交的部分，即

其中:Ainsonif與Abeam的計算公式為

式中:c是聲波在水中的傳播速度，ΘT為發(fā)射陣平行航跡方向波束寬度，ΘR為接收陣垂直航跡方向波束寬度，R為基陣到海底的距離.

系統(tǒng)處理增益包括前置放大器增益、自動增益、時變增益、波束形成處理增益以及接收陣靈敏度級等，可由實際測量獲得.

圖3 有效聲照射面積示意Fig.3 Schematic diagram of effective insonified area

圖4 成像數據處理流程Fig.4 Diagram of imaging data processing

將上述得到的各變量值代入到式(6)便可計算出BS，即地貌信息.但由于反向散射強度有很強的角度關系，尤其是在垂直入射時比較明顯，如果直接成像將無法區(qū)分是海底特征還是入射角度對反向散射強度產生了影響，這會導致成像效果不理想，所以在最終成圖之前需要將其消除，以使反向散射強度只與海底特征有關［10］.由于傳統(tǒng)的Lambert定律模型對垂直入射區(qū)域的角度關系無法合理描述，且對多種海底類型不能準確地預報［11］，所以采用一種經驗修正方法對整個航跡測得的反向散射強度角度趨勢進行平均處理，具體算法表示為式中:BS(θ)為在θ角度處取樣窗內的所有反向散射強度數據，(θ)為上述所有數據的平均值.(β)為角度是β時反向散射強度數據的平均值，加此項的目的是使修正結果可以與其他修正方法更好地相比較，在0°～50°間任何角度都可以作為參考點.通過綜合比較不同角度下現有的試驗數據處理效果，選擇采用40°時的計算值.

將修正后的反向散射強度數據通過量化處理得到各檢測點的灰度值，結合其空間位置信息便可得到海底聲吶圖像.基于上面的分析，利用多波束測深聲吶進行海底成像的數據處理流程如圖4所示.

1.5 算法性能特點

基于多子陣檢測法提出的多波束測深聲吶海底成像技術具有以下性能特點:

1)因為多子陣檢測法可在各波束內得到每個采樣時刻海底回波的DOA估計，所以在此方法基礎上估計的地貌特征能具有較高的分辨率.

2)因為該方法對海底各檢測點的空間位置和回波強度估計是共點測量，所以獲得的地形數據與地貌數據能夠準確地對應，避免了兩者復雜的數據融合問題，進而提高成像質量.

3)因為所得的地形數據量與地貌數據量相同且一一對應，所以該方法便于海底地形地貌信息的四維顯示以及對地物信息的辨識.

2 試驗數據處理

2.1 試驗條件

圖5 多波束測深聲吶系統(tǒng)Fig.5 Diagram of system com position of multibeam bathy metric sonar

為驗證上述方法，利用哈爾濱工程大學自主研制的某型號多波束測深聲吶對吉林省松花湖某區(qū)域采集的湖底回波數據進行處理分析［12］.圖5給出了該多波束測深聲吶系統(tǒng)組成照片，其采用水平放置的“T”型組合聲學基陣，發(fā)射基陣為多元弧陣，接收基陣為由80個陣元組成的均勻線陣，共形成127個接收波束.陣元間距為半波長，發(fā)射信號的工作頻率為300 kHz，采樣頻率為40 kHz，試驗中選用的脈沖寬度為1 ms.

2.2 算法性能驗證

首先，利用傳統(tǒng)的幅度檢測法WMT法與多子陣檢測法對試驗數據進行處理.由于大量的處理結果具有相似的特點，所以這里只取其中一組結果進行說明.圖6(a)、(b)分別是利用這2種方法得到的回波TOA和DOA估計，對比2幅子圖可以看出，WMT法在一個波束內只能獲得一個測量值，而多子陣檢測法可以對海底進行近似連續(xù)的測量，所以在各個波束內可獲得大量的回波強度值在橫向測量區(qū)域密集地分布，從而具有較好的成像分辨能力.前面提到，多子陣檢測法理論上在每一個采樣點上均可獲得海底回波的DOA信息，而實際上由于外側波束的噪聲干擾較大，鏡像波束內的回波相位信息隨采樣時間的變化非常快等原因，在這些區(qū)域會產生一些測量誤差較大的點，因此需相應地將其剔除.

圖7為一組反向散射強度數據的測量值及利用經驗修正方法消除角度關系后的修正結果，從圖中可以看出，在鏡像區(qū)域反向散射強度受到角度的影響較大，而經修正后其角度關系在各角度區(qū)域都能夠較好地得到抑制.

圖6 WMT法與多子陣檢測法對湖試數據的處理結果Fig.6 Results of lake experiment data processed by WMT and the multiple sub-array detection method

圖7 成像數據測量及角度關系修正結果Fig.7 Measurement and correction of the angular dependence of imaging data

2.3 湖底地貌測量

為進一步考核算法的有效性，對松花湖某湖底區(qū)域進行大面積成像測量.圖8為某一航跡成像數據的處理結果，圖中x軸、y軸坐標是將各檢測點的平面大地坐標進行通用橫軸墨托卡(UTM)投影得到［13］，圖像的像素值體現了反向散射強度的大小，顏色越亮代表的強度值越大.

從圖8(a)左下角和右上角的局部放大效果可以看出，該成像測量方法能較細膩地描述湖底地貌的變化.

圖8(b)給出了該航跡湖底地形、地貌信息的一體化表征結果.圖中的z軸坐標是各湖底檢測點相對湖面的深度，很明顯該區(qū)域地形從右上角到左下角起伏變化范圍在50～100 m之間，而地貌信息在該地形圖上得到了一一對應，實現了多波束測深聲吶地形與地貌的共點測量以及一體化表征.

圖8 松花湖某一航跡測區(qū)測量結果Fig.8 Measurement results of a track area in the Songhua Lake

3 結束語

利用由多子陣檢測法獲得的大量海底回波TOA與DOA數據對各海底檢測點的空間位置以及回波強度進行估計，并對其估計結果進行修正.此方法獲得的各檢測點處的反向散射強度數據與其空間位置信息能實現準確地一一對應，避免了2種數據的復雜融合問題，提高了成像質量.湖試數據的處理結果驗證了算法的有效性，并結合修正后數據的成圖效果看出該方法可以獲得較好的空間分辨能力，能較清晰地反映湖底地貌變化，實現地形地貌的共點測量及一體化表征.下一步將利用地貌圖像開展底質分類、海底特征分析等研究工作.

［1］BROWN C J，BLONDEL P.Developments in the application ofmultibeam sonar backscatter for seafloor habitat mapping［J］.Applied Acoustics，2009，70(10):1242-1247.

［2］PRESTON J.Automated acoustic seabed classification of multibeam images of Stanton Banks［J］.Applied Acoustics，2009，70(10):1277-1287.

［3］DEKEYZER R T，BYRNE JS，CASE JD，et al.A comparison of acoustic imagery of sea floor features using a towed side scan sonar and amultibeam echosounder［C］//Oceans＇02 MTS/IEEE.［s.l.］，2002，2:1203-1211.

［4］LINGSCH SC，ROBINSON C S.Acoustic imagery using a multibeam bathymetric system［J］.Marine Geodesy，1992，15(2-3):81-95.

［5］SIMONSD G，SNELLEN M.A Bayesian approach to seafloor classification usingmulti-beam echo-sounder backscatter data［J］.Applied Acoustics，2009，70(10):1258-1268.

［6］LOCKHARTD，SAADE E，WILSON J.New developments inmultibeam backscatter data collection and processing［J］.Marine Technology Society Journal，2001，35(4):46-50.

［7］周天，李海森，朱志德，等.多波束測深系統(tǒng)多子陣檢測法的改進及其性能分析［J］.聲學技術，2005，24(3):152-156.ZHOU Tian，LIHaisen，ZHU Zhide，et al.Performance analysis of improved multiple sub-array detection method in multi-beam bathymetry system［J］.Technical Acoustics，2005，24(3):152-156.

［8］魏玉闊.多波束測深假象消除與動態(tài)空間歸位技術［D］.哈爾濱:哈爾濱工程大學，2011:43-45.WEIYukuo.Technique of bathy metric artifact elimination and seafloor footprint positioning for multibeam bathymetry［D］.Harbin:Harbin Engineering University，2011:43-45.

［9］De OLIVEIRA JA M，HUGHESC JE.Recovering wide angular sector multibeam backscatter to facilitate seafloor classification［C］//United States Hydrographic Conference.Norfolk，VA，2007.

［10］GAVRILOV A N，PARNUM IM.Fluctuations of seafloor backscatter data from multibeam sonar systems［J］.IEEE Journal of Oceanic Engineering，2010，35(2):209-219.

［11］PARNUM IM.Benthic habitat mapping using multibeam sonar systems［D］.Perth:Curtin University of Technology，2008:18-20，101-103.

［12］LIHaisen，YAO Bin，WEIYukuo，et al.Shallow water high resolution multi-beam echo sounder［C］//Oceans 2008-MTS/IEEE.Kobe，2008:1051-1055.

［13］趙建虎，劉經南.多波束測深及圖像數據處理［M］.武漢:武漢大學出版社，2008:124-135.