基于Kongsberg EM 多波束的Snippet数据处理方法

胡 俊吴永亭豆虎林李治远

(自然资源部 第一海洋研究所,山东 青岛266061)

目前,换能器设计、数字电路、信号处理、定位定姿、图形显控等技术的发展使得电子多波束测深系统(Multibeam Echosounder,MBES)的综合性能较20世纪有了飞跃式的提升,多波束已成为海底地形地貌调查的常规水下声学设备[1]。长久以来,随水深一同获取的后向散射回波强度仅仅被当作水深的附属数据,后来研究人员发现后向散射强度与体积非均质性(沉积物粒径、生物扰动、地质分层)和界面粗糙度(海底基质、微地形、地貌等)有关,包含着定性或定量的海底底质类型或组成等信息[2-3],这使得基于多波束的声学底质分类逐渐成为海洋声学底质探测的重要方法和手段[4],并应用于海底生物环境制图[5-6]。多波束后向散射信号是水面、水体、海底以及各种噪声源的综合散射信号,还受多波束系统声源级、发射接收指向性和接收灵敏度影响,为了得到海底的固有散射强度信息,必须消除这些影响。

国内外学者对后向散射强度数据的处理已进行了广泛研究[7-11],但提出的数据处理方法多是针对某一特定项改正,且处理过程中忽视了声呐系统本身的特性,然而对于不同的声呐系统、不同调查目的和不同测量环境下所获取的Snippet数据在处理方法上存在差异,即使是同一艘船上的同一个多波束声呐对同一区域在不同时间和不同航向扫测得到的后向散射强度仍然差别很大。业界缺乏标准化的Snippet数据获取方法和处理流程,且多波束硬件和软件厂商未提供详细的解释文档和资料,这些都限制了后向散射强度数据的进一步利用。为此,本文从声呐方程出发,针对Kongsberg EM 系列多波束系统详细论述了声波从发射到接收过程中经过的一系列处理,给出了一套完整的数据后处理流程,同时对EM 实时改正模型进行了分析和改进,最后利用实测EM2040数据对提出方法进行了验证。

1 Snippet数据处理流程

单位面积海底固有散射强度BS f(β)与3个因素有关:海底底质类型(孔隙度、粗糙度和粒径等)、声波的发射频率f、声波到达海底的入射角β[12];而多波束声呐记录的回波强度RL是原始声呐信号经水体、海底和声呐系统共同作用后的结果,目标散射强度TS不但取决于海底的固有散射强度,还与瞬时反射面积有关[12]:

式中:RL为接收声源级,SL为发射声源级,TL为传播损失,A为声照面积,DT为发射指向性,DR为接收指向性,SH为接收灵敏度,G R为接收增益。因此,多波束后向散射强度处理的本质就是从记录的回波强度RL中恢复单位面积海底固有散射强度BS f(β),这个过程称为辐射校正。

1.1 与声呐系统相关的改正

与声呐系统相关的改正分为接收增益改正和与换能器特性相关的改正(声源级、发射接收指向性和接收灵敏度)。接收增益包括信号AD(数字和模拟)转换之前和之后施加的固定和时变增益(TVG),目的是补偿声波的球面传播损失和水体的吸收损失,以保持回波信号在AD 转换的动态范围内从而实现正确底跟踪。对于不同的多波束系统,其采用的增益方法和实施阶段都不相同,最简单的一种增益方式为[13]

式中:R为斜距;G R的补偿范围比较模糊,主要用于水体图像的实时显示。

多波束采集软件通常会允许用户实时设置参数A,B和C并记录在采集文件中,因此在后处理时需要将Snippet回波信号中的接收增益精确移除,回波信号经过增益移除后的回波强度为

多波束测深仪厂商会提供声源级参数,换能器接收灵敏度、发射和接收指向性可以通过水池实验获得,但会因安装质量、电器元件性能以及使用环境而存在差异,且随着多波束的使用时长而发生改变。通常测量声源级、发射指向性、接收指向性和接收灵敏度的联合作用称为声强校准,得到的结果称为波束指向性模型(Beam Pattern,BP)。单波束和分裂波束鱼探仪可以用反射强度已知的金属(铜或钨)球校准,但是这种方法难以应用于多波束,原因是多波束波束角太小、开角太大、工作频率和模式太多,且校准要求极好的海况。因此目前还没有合适的现场校准方法[14-15],实际作业时仍采用硬件出厂时的BP参数。BP与波束发射接收角和姿态补偿有关:对于无姿态补偿的系统,BP的参考波束角是相对于-3 d B声轴的;而对于有姿态补偿的系统,其是相对于竖直面的。经过改正后的回波强度为

1.2 与水体和海底特性相关的改正

与水体有关的改正TL项分为球面传播损失和水体的吸收损失,均与斜距R有关[13]:

式中:α为吸收系数,其不仅与频率有关,还与温度、压力、盐度、p H 有关;c为声速;t为双程传播时间。

表1为常用声呐频率对应的海水吸收系数。通常多波束采集软件会根据输入的CTD(温盐深)剖面或声速剖面结合平均盐度和发射频率计算出对应深度的平均吸收系数,EM 多波束在不同发射扇区发射频率是变化的,因此不同扇区的吸收系数有差异[16]。

表1 常用声呐频率对应的海水吸收系数[17]Table 1 Absorption coefficients of seawater corresponding to commonly used sonar frequencies[17]

目标散射强度TS是波束脉宽覆盖的区域体散射量,为得到单位面积固有散射强度BS f(β),需进行声照面积改正,其关键在于准确计算出波束打到海底的实际照射面积A[18],经过改正后的回波强度为

EM 多波束为了在达到更大量程的同时保持距离分辨率,在边缘扇区使用了FM(线性调频)信号,FM技术使用匹配滤波的结果会导致脉宽被压缩,其用于声照面积改正的等效脉宽τeff为

式中:BW为对应扇区的带宽,并且需要考虑由于匹配滤波而增加的信噪比:

式中:T为对应扇区的脉冲发射时间,通常为第一个脉冲的发射时间加上扇区时延。因此实际测量时为了获得高质量的图像数据尽量使用CW 信号代替FM 信号。此外,还需要考虑地形的影响,Lurton用实验验证了在局部微观地形的作用下,从不同方向扫测得到的后向散射强度差异最大可达到11 dB[19];然而基于相邻波束点坐标计算得到的地形倾角误差太大,通过导入水深DEM 估计航迹向地形倾角和垂直航迹向地形倾角是Caris等后处理软件所采用的方法。本文在此基础上进行了改进、直接对不经过地形改正的RL3回波信号采用最大脉宽τ对应的分辨率cτ/2进行基于扫描填充法的地理编码格网成图,按照距离倒数加权动态更新每个格网的水深信息和强度信息,待所有Snippet点扫描完毕,计算每个格网点的坡度信息并进行地形倾角改正。

1.3 角度响应改正

Snippet数据处理结果分为回波强度-角度响应曲线和去角度响应的地理编码图像,前者主要用于ARA(Angular Range Analysis)反演底质类型,而地理编码图像对灰度均一性有较高的要求。经过辐射校正后的回波强度表现出与波束入射角相关的特点,比较常见的现象就是中央波束在图像上呈现明显的亮区异常,因此将回波强度表达为与入射角无关的量的过程称为角度响应改正。实际上角度响应特征是不同海底底质类型的固有特性,无法通过模型完全消除,比较合适的做法是将不同入射角下的回波强度归一化到某一参考角度或某一参考角度范围内的平均回波强度,即AVG 改正,具体改正步骤:

1)纵向选择连续的底质类型一致的子区数据段;

2)横向上基于一定的角度分辨率(通常为1°)统计每个角度下对应的平均回波强度得到角度-强度查找表;

3)在角度-强度查找表中选择某一参考角度或某一参考角度范围内的平均强度作为整个区段的参考回波强度基准值;

4)将原始回波强度减去角度-强度查找表中的平均强度再加上整个区段的参考回波强度基准值得到最终的回波强度。

其中步骤2)得到的角度-强度查找表也可以用拟合的经验模型曲线代替[20]或分段拟合直线代替[21]。此方法的难点在于分段长度的选取:选择较长的连续ping数据段能够降低统计误差,增加横向上的平滑性,但是容易将不同的底质类型混合造成图像纵向上的“假象”;而选择较短的ping数据段则会导致横向上灰度不连续。通常做法是统计改正ping前后各一段的角度响应曲线进行移动平均。声照面积改正、地形倾角改正、BP 改正均与入射角相关,其改正残余量会同BS f(β)耦合在一起共同参与AVG 改正,但是与入射角无关的不正确的固定增益和过度的TVG 增益会导致AVG 改正不彻底。与扫描发射多波束(如Sea Beam 系列多波束)不同的是EM 多波束采用扇区发射的方式导致其不能直接采用基于全波束角范围内的AVG 改正,需要考虑相邻扇区之间的声源级差异,因此实际改正时需要分扇区进行角度响应统计和改正。

将Snippet数据处理流程分为4步3阶段(图1):通过对原始数据解码获得初始回波强度RL结合声呐系统特性进行相关改正得到与声呐系统无关的回波强度RL1;再通过物理模型改正去掉水体和海底特性的影响得到与入射角相关的回波强度RL2;利用RL2进行角度响应分析和改正得到某一参考角度下的回波强度RL3,最后地理编码成图并去掉地形的影响。

图1 Snippet数据处理流程Fig.1 Flow chart of snippet data processing

1.4 EM 模型改进

对于Kongsberg EM 多波束,all文件记录的Snippet回波强度已经过声源级、接收灵敏度、接收增益、发射接收指向性、水体吸收损失、声照面积和AVG 改正。EM 多波束在采集端使用的实时改正模型[18]:

式中:BS N为垂直入射时的平均后向散射强度,BS0为波束角等于T(镜面反射和Lambert反射临界角)时的平均后向散射强度,其中k=1/cosT,s=R/R I,R I为波束垂直入射时的距离,R为斜入射距离,T可在EM多波束采集软件SIS中手动设置,默认值为6°。此外,EM 多波束使用了内置的BP改正模型,相关参数记录在甲板单元PU 的BSCorr.txt文件中(可通过ftp下载或修改该文件)。

EM 多波束在采集端进行了所有改正,其得到的后向散射强度近似接近于海底固有散射强度。因此,通常情况下EM 多波束用户可以直接利用all文件的图像数据而不经过任何改正;但是EM 的改正模型过于简单,难以适用于复杂的底质类型和特殊地形,因此进行了改进:

1)移除接收增益GO1,该参数记录在all文件的安装包中;

2)根据底检测质量因子,对底跟踪失败导致斜距R估算错误的波束依据式(5)重新做传播损失和吸收损失修正;

3)EM 声照面积改正忽略了地形的影响,按照本文1.2节中提出的方法增加地形倾角改正;

4)AVG 改正中,EM 将整个反射区分为小入射角区(线性模型)和漫反射区(Lambert模型),由于是实时改正,EM 只能基于前向数据估计垂直入射和斜入射的参考回波强度,因此需要移除EM 模型中的Lambert改正项,然后基于前后ping分扇区按照1.3节中的方法进行角度响应改正。

至于BP改正,若无条件做声强校准(Amplitude Calibration)得到准确的BP模型,不建议重新改正。

2 试验与分析

为验证本文提出方法的实用性,选取Kongsberg浅水多波束EM2040采集的all格式数据应用本文提出的方法进行处理;EM2040数据选自Kongsberg官方示例的沉船数据和黄海的一个被打捞过的沉船区域数据。EM 官方示例数据的沉船位于测线航迹正下方,整个条带水深范围为23.8～25.2 m,地形较为平坦,沉船的挤压导致附近区域地形略微下陷。所用2040波束脉宽为50μs,Lambert临界角为10°,经过解码、水深融合等一系列常规处理后,采用0.055 m 的分辨率直接对原始数据进行扫描填充的地理编码法生成图像(图2a),图2a中沉船清晰可见。但是EM 不准确的实时改正模型造成了图2a和2b中沿航迹方向在左右舷10°附近出现连续的“暗线”异常,沉船的回波强度(区段2)与地形的回波强度(区段1和3)有明显的差异,从图3中的AR(Angular Response)曲线可以看到区段1和3的回波强度保持在-30 dB上下的水平,而区段2的回波强度范围为-45～-10 dB,导致在船尾和船首附近区域横向上的灰度过度不均衡,图像上呈现出明显的分界线。为此,通过移除EM 的接收增益和lambert模型,利用地理编码时生成的等分辨率DEM 做地形改正及传播损失和吸收损失修正,然后重新成图得到图2b,实现了整个图像横向上的灰度均衡。但是天底方向的“暗线”仍然存在,采用前后各15 ping的移动窗分3个扇区对图像做AVG 改正(图2c),由图2c可知图2a和图2b中的“暗线”得到消除,而且沉船附近的灰度与整个条带灰度更加一致。

图2 EM2040官方沉船数据Snippet图像处理结果对比Fig.2 Comparison of processing results of Snippet images of official sunken ship data collected by EM2040

为了说明改正效果,统计了区段1和区段2改正前后的角度响应曲线,如图3b和3c所示,原始回波强度在左右舷10°和45°(边缘扇区和中央扇区分界角)附近存在“跳变”,经过分扇区AVG 改正后该现象得以消除,回波强度在整个开角范围内趋于一致。表2为原始数据成图和校正后图像的平均回波强度、标准差和熵值对比,校正后的图像保持了原来的回波强度水平,但是标准差明显减小,整个图像的回波强度分布更加均匀。而熵值的减小是由于AVG 改正的平均效应模糊了图像细节,这表明AVG 改正后的图像不利于目标识别,但是能显著消除或削弱入射角对回波强度的影响。

表2 原始数据所成图像和校正后图像对比Table 2 Comparison of images generated from original data and corrected image

图3 原始数据和校正后数据的角度响应曲线Fig.3 Angular response curves of original data and corrected data

为验证本文所提出的方法在多条带中的应用,将黄海某一区域的EM2040校准数据按照图1 中的方法进行处理得到图4,由图4中可知,本文提出的方法不仅能有效地消除中央区域“亮区”异常和扇区间的灰度不均衡,还能削弱拼接痕迹。在处理过程中发现:AVG 改正时波束入射角需顾及横摇姿态角的影响,否则会导致图像边缘部分产生畸变;地理编码时可同时生成等分辨率的DEM,视地形平坦程度决定是否做地形校正。Kongsberg 目前还未给用户提供EM2040的内置BP 模型,因此BP的修正只能与AVG 改正一同处理。

图4 多条带Snippet图像处理结果对比Fig.4 Comparison of processing results from multiple swaths

3 结 语

多波束后向散射强度的处理不仅需要顾及声呐系统本身的特性,还需要考虑到水体介质的作用以及海底目标的散射特征。不同多波束硬件系统有各自的后向散射强度处理流程和方法,本文以Kongsberg EM多波束为例,提出了一套标准化的Snippet数据处理流程,将数据处理分为采集、后处理、应用三个层次,将回波强度的处理水平分为4个阶段。针对EM 实时改正模型无法应对复杂的海底底质类型和特殊地形的情形,提出了相应的后处理改进方法和建议,并对典型EM2040实测数据进行了处理,结果表明本文的方法能够有效消除Snippet回波信号中的辐射畸变,显著提高图像质量。受多波束系统BP改正不完善和AVG 校正不彻底的限制,目前还没有更好的办法能直接获取反应海底底质类型的纯散射强度,所生成的图像都是基于某一特定入射角或某一入射角区间的参考回波强度。尽管我们获得了大量的海底散射数据,但是却无法统一量化标准,只能寄希望于未来的多波束硬件和软件在处理流程和方法上进行改进,并找到合适的多波束声强校准方法。

参考文献(References):

[1]PRATSON L F,EDWARDS M H.Introduction to advances in seafloor mapping using sidescan sonar and multibeam bathymetry data[J].Marine Geophysical Researches,1996,18(6):601-605.

[2]CLARKE J E H,MAYER L A,WELLS D E.Shallow-water imaging multibeam sonars:a new tool for investigating seafloor processes in the coastal zone and on the continental shelf[J].Marine Geophysical Researches,1996,18(6):607-629.

[3]CLARKE J E H,DANFORTH B W,VALENTINE P.Areal seabed classification using backscatter angular response at 95 k Hz[C]∥NATO SACLANTCEN.Conference on High Frequency Acoustics in Shallow Water.Italy:Lerici,1997:243-250.

[4]TANG Q H,JI X,DING J S,et al.Research progress and prospect of acoustic seabed classification using multibeam echo sounder[J].Advances in Marine Science,2019,37(1):1-10.唐秋华,纪雪,丁继胜,等.多波束声学底质分类研究进展与展望[J].海洋科学进展,2019,37(1):1-10.

[5]LAMARCHE G,ORPIN A R,MITCHELL J S,et al.Benthic habitat mapping[M]∥Hoboken,U.S.:John Wiley&Sons,Ltd,2016,80-102.

[6]HASAN R C,IERODIACONOU D,LAURENSON L.Combining angular response classification and backscatter imagery segmentation for benthic biological habitat mapping[J].Estuarine Coastal&Shelf Science,2012,97:1-9.

[7]JIN S H,ZHAI J S,LIU Y C,et al.Data refinement of Simrad EM multibeam backscattering intensity[J].Science of surveying and mapping,2010,35(2):106-108.金绍华,翟京生,刘雁春,等.Simrad EM 多波束反向散射强度数据精处理研究[J].测绘科学,2010,35(2):106-108.

[8]TANG Q H,ZHOU X H,DING J S,et al.Data processing of multibeam backscattering intensity[J].Haiyang Xuebao,2006,28(2):51-55.唐秋华,周兴华,丁继胜,等.多波束反向散射强度数据处理研究[J].海洋学报,2006,28(2):51-55.

[9]JIN S H,ZHAI J S,LIU Y C,et al.Influence of seafloor incidence angle on multibeam backscatter intensity and corrected method[J].Geomatics and Information Science of Wuhan University,2011,36(9):1081-1084.金绍华,翟京生,刘雁春,等.海底入射角对多波束反向散射强度的影响及其改正[J].武汉大学学报(信息科学版),2011,36(9):1081-1084.

[10]HELLEQUIN L,LURTON X,AUGUSTIN J M.Postprocessing and signal corrections for multibeam echo sounder images[C]∥OCEANS′97 MTS/IEEE Conference Proceedings,1997(21):23-26.

[11]AUGUSTIN J M,LULEBLANC E,FOSTER B.New algorithm for multibeam imagery processing[J].Marine technology reporter white papers,hydrographic edition,2015:40-44.

[12]AUGUSTIN J M,LURTON X.Image amplitude calibration and processing for seafloor mapping sonars[C]∥Oceans.Europe,20-23 June 2005(1):698-701.

[13]Kongsberg Maritime AS.Kongsberg EM Series Multibeam Echo Sounder Datagram formats[EB/OL].Norway:Kongsberg Simrad Company,2015(2017-01-02)[2019-08-10].https://www.kongsberg.com/globalassets/maritime/km-products/product-documents/160692_em_datagram_formats.pdf.

[14]HEATON J L,RICE G,WEBER T C.An extended surface target for high-frequency multibeam echo sounder calibration[J].The Journal of the Acoustical Society of America,2017,141(4):388-394.

[15]LANZONI C,WEBER T.A method for field calibration of a multibeam echo sounder[J].Oceans,2011:1-7.

[16]CARVALHO R,CLARKE J E H.Proper environmental reduction for attenuation in multi sector sonars[C]∥2013 IEEE/OES Acoustics in Underwater Geosciences Symposium.Rio de Janeiro,Brazil:IEEE,2015:1-6.

[17]AINSLIE M A,MCCOLM J G.A simplified formula for viscous and chemical absorption in sea water[J].Journal of the Acoustic Society of America,1998,103(3):1671-1672.

[18]HAMMERSTAD E.Backscattering and seabed image reflectivity[EB/OL].Norway:Kongsberg Simrad Company,2000[2019-08-10].https://www.kongsberg.com/globalassets/maritime/km-products/product-documents/backscattering-and-seabed-image-reflectivity.

[19]LURTON X,ELEFTHERAKIS D,AUGUSTIN J M.Analysis of seafloor backscatter strength dependence on the survey azimuth using multibeam echosounder data[J].Marine Geophysical Research,2018,39:183-203.

[20]JIN S H,XIAO F M,BIAN G,et al.A method for extracting seabed feature parameters based on the angular response curve of multibeam backscatter strength[J].Geomatics and Information Science of Wuhan University,2014,39(12):1493-1498.金绍华,肖付民,边刚,等.利用多波束反向散射强度角度响应曲线的底质特征参数提取算法[J].武汉大学学报(信息科学版),2014,39(12):1493-1498.

[21]YAN J,ZHANG H M,ZHAO J H,et al.Improved algorithm for angular response model of multi-beam sonar backscatter data[J].Journal of Surveying and Mapping,2016(11):1301-1307.严俊,张红梅,赵建虎,等.多波束声呐后向散射数据角度响应模型的改进算法[J].测绘学报,2016(11):1301-1307.