文章信息
- 付桂合, 马鑫程, 冯成凯, 印萍, 梅赛, 阳凡林. 2020.
- FU Gui-he, MA Xin-cheng, FENG Cheng-kai, YIN Ping, MEI Sai, YANG Fan-lin. 2020.
- 混浊水域多波束反向散射强度传播损失改正
- Transmission-loss correction of multibeam sonar backscatter strength in turbid waters
- 海洋科学, 44(5): 107-114
- Marina Sciences, 44(5): 107-114.
- http://dx.doi.org/10.11759/hykx20190424002
-
文章历史
- 收稿日期:2019-04-24
- 修回日期:2019-10-28
2. 中交公路规划设计院有限公司, 北京 100088;
3. 青岛海洋地质研究所, 山东 青岛 266071;
4. 自然资源部海岛(礁)测绘技术重点实验室, 山东 青岛 266590
2. CCCC Highway Consultants CO. Ltd, Beijing 100088, China;
3. Qingdao Institute of Marine Geology, Qingdao 266071, China;
4. Key Laboratory of Surveying and Mapping Technology on Island and Reef, Ministry of Natural Resources, Qingdao 266590, China
多波束声呐(multibeam sonar)是目前使用最为广泛的海洋测量仪器之一, 不仅能够获得水深数据, 还能够获得反向散射强度(backscatter strength, BS)。反向散射强度经过一系列的处理后, 可以绘制海底声呐图像, 为海底底质分类、海洋目标探测、海洋环境调查等诸多领域提供基础数据支撑[1-2]。多波束反向散射强度用于定量描述目标反射特性, 不仅受海底底质的影响, 而且还依赖于仪器特性及海洋环境[3-4]。多波束回波数据处理过程中, 去除海洋环境影响尤为重要, 海洋中的诸多因素都对声信号的传播产生影响, 这些因素主要包括声速剖面、水体温度、水体盐度、压强以及悬浮颗粒物等[5-6]。
国内外学者[7]都已展开多波束回波声传播损失的研究, 现有的反向散射强度传播损失改正方法主要有多波束换能器增益处理和声学模型后处理, 多波束换能器增益是通过时间变化增益(time varied gain, TVG)减少传播损失的影响, 唐秋华等[8]分析了强度数据获取过程中系统进行的实时补偿。Beaudoi等[9]分析了多波束系统换能器接收增益对回波强度值的影响, 声学模型法主要是通过声传播公式推导回波强度与海洋环境参数的关系, 彭临慧等[5]研究了混浊海水对声传播的影响, Alexandre等[4]对多波束反向散射强度传播损失与海洋环境参数的关系进行了相关研究。
在近岸河海交汇处, 河流携带大量的泥沙汇入海洋, 产生大片混浊海域。混浊海水中的悬浮颗粒物, 通过热传导以及黏滞作用导致声波衰减, 对多波束声呐系统的工作性能产生显著影响, 降低了声呐的识别能力。因此研究混浊水域中的多波束传播损失成为关键问题。波束传播损失与波束的传播路径密切相关, 未考虑声速剖面对波束传播路径的影响。而传统的传播损失与波束的传播路径密切相关, TVG改正方法忽略了海水介质中环境参数的动态性, 而传统的多波束传播损失声学模型通常只考虑了纯净海水黏滞性对波束传播的影响[10], 未考虑混浊海水中悬浮颗粒物对海水黏滞性的影响, 在复杂海洋环境下适应性较弱, 改正效果不理想。
针对上述问题, 本文从声学机理出发, 考虑混浊水域中悬浮颗粒物对黏滞吸收的作用, 提出了混浊水域多波束声传播损失改正模型, 针对海洋环境参数的动态性, 对每个波束分层计算传播损失, 最后将各层传播损失累计, 获得每个波束的传播损失。
1 多波束传播损失影响因素 1.1 声传播基本原理波束发射后, 经海水传播, 再到接收的整个过程形成了声呐方程。声呐方程的影响因素主要包括声传播介质、海底及水体目标、背景干扰和声呐参数。图 1展示了波束信号从发射到接收的过程, 图 1中各个参数表示如下: SL为多波束系统的声源级; TL为波束的传播损失; TS为目标强度; A为声照区面积; EL为回波能级; NL为海洋环境的声噪水平; A/D为模数变换; Dt和Dr分别发射和接受指向性; 为Gb为ADC阶段之前实施的增益; Ga为ADC阶段之后实施的增益; RL为系统增益之后接收回波能级。BSB为海底反向散射强度。则声呐方程可表示为[10]:
$ \left\{ \begin{align} & \text{EL=SL}-\text{2TL+TS}-\text{NL+}{{\text{D}}_{\text{r}}} \\ & \text{TS=B}{{\text{S}}_{\text{B}}}\text{+10}\lg A \\ \end{align} \right., $ | (1) |
因此根据声呐方程, 影响海底固有反向散射强度BSB的因素主要有SL、EL、TL、NL、Dr等。传播损失TL定量描述了声波传播一定距离后声强度的衰减变化, 传播损失由扩展损失和吸收损失两部分构成, 图 2展示了影响传播损失的主要因素。扩展损失是由声波本身的传播特性决定的, 声波在传播过程中波阵面不断扩大, 单位面积能量减少, 衰减的速率与波阵面的表面积成正比。对于多波束声呐系统发射的球面波, 其衰减的速率与传播距离平方成正比[11]。吸收损失的本质是声传播过程中部分声能转化为其他能量(如热能), 且不可逆, 主要包括化学弛豫吸收和黏滞吸收, 弛豫吸收是指海水中的分子在声波作用下发生离解和缔合, 消耗声波的能量, 黏滞吸收是由于惯性黏滞作用将速度梯度能量转化为热能而产生的能量损失[12]。
1.2 常用传播损失改正方法 1.2.1 TVG改正在水深测量过程中多波束系统会自动添加固定增益(Fixed Gain)及时间变化增益, 即根据信号发射后的时间长短对返回信号进行放大处理[7]:
$ {{\text{T}}_{tvg}}\text{ }=\text{ }{{S}_{\text{p}}}\lg R+2\alpha R\text{/1 000}+G, $ | (2) |
式中, α为吸收系数(dB/km), R为斜距(m), Sp为扩展损失系数, G为固定增益。
时间增益是另一种形式的传播损失改正, 吸收系数α和传播损失系数Sp是用户在进行初始设置时输入的参数值, 与测区真实的吸收系数是不同的, 存在着一定的不确定性和局限性, 因此需要对采集到的回波强度进行去除TVG改正, 在数据后处理中重新计算传播损失。
1.2.2 声传播损失计算模型传统的多波束声传播损失由两部分组成, 一部分是由波束波阵面随距离扩展产生的声强衰减(扩展损失), 一部分是由水的黏滞性和水中二价盐等化学弛豫过程引起的声强衰减(吸收损失)[10]:
$ \text{TL}=20\lg R+\alpha R/1\ 000, $ | (3) |
式中, R为射程(m), α为吸收系数(dB/km)。α是声呐频率f和传播介质特性(盐度、温度、深度和pH值)的函数[13]:
$ a=\frac{{{A}_{\text{1}}}{{f}_{1}}f}{{{f}^{2}}+f_{1}^{2}}+\frac{{{A}_{\text{2}}}{{p}_{2}}f_{2}^{2}{{f}^{2}}}{{{f}^{2}}+f_{2}^{2}}+{{A}_{\text{3}}}{{p}_{3}}{{f}^{2}}, $ | (4) |
式中, f频率(kHz)、T为温度(℃)、S为盐度和Z为水深(km), c为声速(m/s), pH值为7.6~8.2级, A1和A2两项表示由硼酸和硫酸镁引起的弛豫衰减。A3项表示由纯水中的黏性引起的衰减。
2 混浊水域声传播损失模型与改正在进行多波束反向散射数据处理时, 发现根据传统方法处理的结果中反向散射强度值随传播距离的增大而减少, 这是由于传统的改正方法未考虑混浊海水中颗粒物对传播损失的影响, 没有完全去除传播损失。因此, 如何处理悬浮颗粒物对吸收损失的影响是混浊海水传播损失计算的关键。本文从声学机理出发, 重新构建了多波束声呐数据吸收损失计算模型, 使用水体温度对海洋环境参数进行自适应分层, 然后基于声速, 对每个波束分层计算吸收系数和传播距离, 最后将各层传播损失累计, 得到每个波束的传播损失, 改正流程图如图 3所示。
2.1 声波传播损失模型改正相较于清澈海水, 混浊海水含有大量的悬浮颗粒物, 影响海水的黏滞性。为此, 本文算法中根据海水混浊程度重新计算黏性引起的衰减αv [14]:
$ \left\{ \begin{align} & {{\alpha }_{v}}= (10\lg {{\text{e}}^{2}})\frac{\varepsilon k{{(\sigma -1)}^{2}}}{2}\frac{\tau }{{{\tau }^{2}}+{{(\sigma +\delta )}^{2}}} \\ & \sigma =\frac{1}{2}+\frac{9}{4\beta \alpha }, \tau =\frac{9}{2\beta \alpha }\left[ 1+\frac{1}{\beta \alpha } \right] \\ \end{align} \right., $ | (5) |
式中, ε=M/ρ′是粒子的体积分数, M为悬浮颗粒物的质量浓度; k=ω/c是波数; c是海水中声速; δ=ρ′/ρ是泥沙颗粒密度与液体密度之比。α是粒子的半径;
从式(5)中可以看出, 黏滞声吸收无法通过水体温度、盐度、深度直接获得, 但可以通过海水动黏滞率与水体温度、盐度、深度建立关系[15]:
$ v=\frac{0.1\left( \sum\limits_{i}{{{p}^{i}}\sum\limits_{j}{{{Q}_{ij}}{{t}^{j}}}}+S\sum\limits_{k}{{{R}_{k}}{{t}^{k}}} \right)}{\rho }. $ | (6) |
考虑海水中悬浮颗粒物造成的黏滞损失, 综合考虑式(4)中的扩展损失、弛豫损失, 则改进的传播损失声学模型为:
$ \text{TL}=40\lg R+2{{\alpha }_{w}}R+2{{\alpha }_{v}}R,$ | (7) |
$ {{a}_{w}}=\frac{{{A}_{\text{1}}}{{f}_{1}}f}{{{f}^{2}}+f_{1}^{2}}+\frac{{{A}_{\text{2}}}{{p}_{2}}f_{2}^{2}{{f}^{2}}}{{{f}^{2}}+f_{2}^{2}}. $ | (8) |
式(7)中, 第一项为扩展损失, 第二项为海水中分子弛豫衰减损失, 第三项为海水黏滞声吸收损失。
2.2 分层计算传播损失不同深度的海洋环境参数常不相同, 而这些海洋环境参数对多波束反向散射强度传播过程的吸收损失有着直接的关系, 如果使用统一的系数进行计算, 会增加误差。因此, 需根据不同深度的海水参数, 对水体进行分层, 计算声线路径上变化的吸收系数, 依次计算每层的传播损失, 直到到达水底, 获得累计吸收损失。首先根据水体温度剖面对波束路径进行分层, 然后利用各层的海洋参数计算传播距离与吸收系数, 计算每一层的吸收损失与扩展损失, 最后获得累计的传播损失。
温盐深培面仪(Conductivity-Temperature-Depth profiler, TCD)或声速剖面仪可获取不同深度海洋环境参数, 现场实测的数据往往采样间隔较小, 获取数据较多, 就浅水多波束系统而言, 多波束每秒能够发射几十次声脉冲(Ping)数据, 每Ping数据中包含几百个波束, 如果采样间隔过小, 算法运算量较大, 对计算效率提出了很高要求, 因此, 需经过数据合理筛选才能使用。仅凭经验等距离进行分段划分, 缺乏分层合理性的判断依据。本文在道格拉斯(Douglas-Peucke)算法[16]基础上, 提出一种适用于多波束传播损失改正的自适应分层方法。
图 4表示的是2016年4月于浙江舟山海域部分实测温度和盐度剖面数据, 水深在30 m左右, 由图中可以看出, 盐度的变换范围在0.5, 对声速影响较小, 因此, 本文算法只考虑使用水体温度进行自适应分层, 具体过程如下:
(1) 利用仪器采集到的数据往往是按等深度间隔分布, 可将离散温度剖面数据表示为(Ti, Di), i = 1, 2, …N。
(2) 如图 5所示,选取温度剖面第一个点(T1, D1)和最后一个点(Tn, Dn), 将两端点连成一条直线, 计算声速剖面曲线上各离散声剖点(Ti, Di)到该直线的距离Zi, 设定阈值为δ, 若{Zi}max < δ, 则保留两端点舍去所有中间点, 若{Zi}max≥δ, 则保留温度剖面中到直线距离最大的中间点, 并以该点为基准将温度剖面分为两部分。
(3) 对上述两部分数据重复过程(2), 直到没有点被舍去为止, 最终将结果分为自适应分层的节点。
声速剖面是水深测量不可忽略的环境变量, 在进行声速测量时, 一些类型的仪器, 直接利用水柱测量声速, 只记录声速和对应的深度, 无法直接获得盐度和温度数据。对于水深测量, 此种类型的数据是有效的。但是对于反向散射强度测量, 需要使用温度和盐度计算吸收系数剖面。
为了克服缺乏温度和盐度信息数据的问题, 在上述内容中提到, 在进行浅水区域测量时, 盐度变化较小, 水体中声速随温度、盐度以及深度的变化而变化, 利用这种相关性可以建立声速经验模型, 在浅水区域, 可以通过式(9)间接确定声速[17]。假设盐度在空间或调查时间内没有显著变化, 使用一定区域内的平均盐度, 根据声速剖面数据和声速经验模型, 得到不同深度的水体温度, 随后使用水体温度进行自适应分层。为验证声速构建吸收系数剖面的精度, 分别使用CTD水体温度数据和声速剖面数据构建了吸收系数剖面, 从图 6中可以看出, 声速和CTD数据获得吸收系数剖面在浅水区域整体趋势相同, 对传播损失的改正影响较小。
多波束测量中, 波束发射存在着一定的角度, 波束通常为非垂直入射, 不同海洋环境下的传播路径完全不同。如果认为声线在水柱中按直线传播, 使用表层声速和旅行时间计算传播距离, 则无法真实还原声速传播路径。
$ \begin{array}{l} c = 1448.96 + 4.591T - 5.304 \times {10^{ - 2}}{T^2} + 2.374 \times \\ {10^{ - 4}}{T^3} + 1.340(s - 35) + 1.630 \times {10^{ - 2}}D + 1.675 \times \\ {10^{ - 7}}{D^2} - 1.025 \times {10^{ - 2}}T\left( {s - 35} \right) - 7.139 \times {10^{ - 13}}{\rm{T}}{{\rm{D}}^3} \end{array} $ | (9) |
本文在常梯度声线跟踪算法的基础上[10], 使用自适应分层的节点, 沿波束传播路径逐层计算传播距离与传播损失。具体每个波束传播损失计算如下:
(1) 根据不同深度的海洋环境数据获得吸收系数剖面;
(2) 从换能器表面开始追加水层, 计算每个层的传播时间和传播距离;
(3) 累计各层的传播时间和传播损失, 累计时间与波束的实测时间比较, 并判断是否追加水层。若时间相同, 则累计的传播损失即波束传播损失, 若累计时间小于实测时间, 回到步骤(2)。若累计时间大于实测时间, 则多追加了传播损失, 需计算多余的传播损失:
$ \Delta {\rm{TL'}} = \alpha 'c'\Delta t',$ | (10) |
式中, α′为层内吸收系数, c′为水层下界声速, Δt′为累计时间与实测时间差值。当测量深度大于获得的声速剖面深度时, 没有海洋环境参数数据, 则利用最后一层数据计算多余部分。
3 实验与分析为验证本文算法, 选取2016年4月于浙江某海域部分实测多波束反向散射数据进行分析。采用R2Sonic 2024多波束测深系统, 其波束开角为120°, 声呐频率为400 kHz, 波束个数为256个, 实验区水深为20~30 m, 位于长江入海口, 大量泥沙流入海中, 故海水混浊。
为减少深度对强度的影响, 实验选取深度较为均匀的3条相邻测线, 利用图 3所示算法流程进行了传播损失改正。为便于比较, 选取未经传播损失改正、TVG改正及传统声学模型改正与本文算法处理结果进行对比分析。在实际处理中, 由于原始反向散射强度数据存在随机误差, 故对连续Ping的反向散射强度取平均, 对数据进行平滑, 消除随机误差影响。将获得CTD数据按照水体温度进行自适应分层, 设定阈值为0.2, 3条测线总共有370万个波束, 对每个波束分层改正传播损失。对原始数据进行声呐参数、声照区面积、入射角效应改正, 保证剩余数据主要受传播损失的影响。
3.1 实验结果对处理后的3条测线数据进行地理编码, 得到测区海底多波束声呐图像(图 7)。从图 7a可以看出, 传播损失对回波强度有较大影响, 声呐图像呈现整体的不均衡性, 很难直接进行底质分类。TVG改正后, 一定程度上去除了传播损失, 但中央区域存在“亮线”, 仍有明显的拼接痕迹; 从图 7c中可以看出, 经传统模型改正后, 整体声呐图像较为均衡。为了更好地对比, 提取图 7中红色方框区域数据, 得到传统模型和本文算法改正局部区域效果图(图 8), 从图 8中可以看出, 经传统模型改正后, 仍然可以看到测线边缘的痕迹; 经本文算法处理后, 有效去除了传播损失, 声呐图像清晰, 测线过渡自然。
3.2 实验分析
由声呐方程可以看出, 回波强度EL是目标强度与传播损失的组合, 若传播损失改正不完全, 强度数据会与传播距离有一定的相关性。Spearman等级相关系数[18]可以评价两个统计变量的相关性, 可以很好地检测变量间变化趋势的相关性:
$ \theta = \frac{{\sum {\left( {{R_i} - \bar R} \right)} \left( {{S_i} - \bar S} \right)}}{{\sqrt {\sum {{{\left( {{R_i} - \bar R} \right)}^2}} {{\left( {{S_i} - \bar S} \right)}^2}} }}. $ | (11) |
Spearman等级相关系数值越接近0, 说明反向散射强度与传播距离的相关性越小, 进而说明传播损失效果越好。提取上述测区中部分反向散射数据, 具体位置为图 7中红色方框所选区域, 计算出反向散射强度和传播距离的Spearman等级相关系数绝对值, 数据改正前、TVG改正、传统模型改正和本文算法改正的Spearman相关系数分别为0.93、0.74、0.44和0.04。通过比较发现, 本文方法Spearman等级相关系数绝对值最小, 说明本文方法传播损失改正效果最佳。
为直观表达, 绘制出强度数据散点图, 并拟合其趋势线, 结果如图 9所示。从图中可以看出, 3种方法均在一定程度上去除了传播损失。由于TVG改正固定声波吸收系数, 反向散射强度仍与传播距离有很大的相关性; 相较于TVG改正, 传统声传播损失计算模型对于传播损失的去除效果较好, 但由于未考虑悬浮颗粒物对传播损失的影响, 随传播距离的增大, 仍具有一定的误差; 本文算法趋势线基本呈直线, 传播损失去除效果最好。
4 结论本文综合分析了混浊水域对多波束声波传播的影响, 提出了混浊海水传播损失改正的方法。通过分析混浊水域悬浮颗粒物对传播损失的影响, 对传统的声学模型进行了改进。并通过实验分析验证了该方法的有效性, 可得到以下结论:
1) 通过对实测数据进行处理分析, 有效去除了多波束回波强度的传播损失, 解决了其造成声呐图像整体明暗不一的问题。经本文算法处理后, 反向散射强度与传播距离拟合趋势线基本呈直线, 反向散射强度与传播距离的相关性明显降低, 海底声呐图像清晰, 满足工程实际需要。
2) 由于声学模型改正方法需要不同深度的海洋环境参数或声速, 并且需要分层计算, 若不能有效减少分层, 计算量较大。本文算法根据水体温度数据进行自适应分层, 在保证精度的情况下, 提高了运算效率。
3) 基于声速, 准确还原波束的传播路径, 并沿波束传播路径分层计算传播距离与传播损失, 大大提高了传播损失的计算精度。
[1] |
唐秋华, 纪雪, 丁继胜, 等. 多波束声学底质分类研究进展与展望[J]. 海洋科学进展, 2019, 37(1): 1-10. Tang Qiuhua, Ji Xue, Ding Jisheng, et al. Research progress and prospect of acoustic seabed classification using multibeam echo sounder[J]. Advances in Marine Science, 2019, 37(1): 1-10. |
[2] |
吴自银, 阳凡林, 李守军, 等. 高分辨率海底地形地貌:可视计算与科学应用[M]. 北京: 科学出版社, 2017: 221-228. Wu Ziyin, Yang Fanlin, Li Shoujun, et al. High-Resolution Submarine Geomorphology:Visual Computation and Scientific Applications[M]. Beijing: Science Press, 2017: 221-228. |
[3] |
Preston J, Christney A, Bloomer S, et al. Seabed classification of multibeam sonar images[C]//IEEE. MTS/IEEE Conference and Exhibition (4). Honolulu: IEEE, 2001: 2616-2623.
|
[4] |
Schimel A, Beaudoin J, Parnum I, et al. Multibeam sonar backscatter data processing[J]. Marine Geophysical Research, 2018, 39(1-2): 121-137. |
[5] |
彭临慧, 王桂波. 中国近海悬浮颗粒物海水声波衰减[J]. 声学学报, 2008, 33(5): 389-395. Peng Linhui, Wang Guibo. Sound attenuation in suspended particulate matter seawater of Chinese sea offshore[J]. Acta Acustica, 2008, 33(5): 389-395. |
[6] |
Brown N, Leighton T, Richards S, et al. Measurement of viscous sound absorption at 50-150 kHz in a model turbid environment[J]. The Journal of the Acoustical Society of America, 1998, 104(4): 2114-2120. |
[7] |
Foote K. Standard-target calibration of active sonars used to measure scattering:principles and illustrative protocols[J]. IEEE Journal of Oceanic Engineering, 2018, 43(3): 749-763. |
[8] |
唐秋华, 周兴华, 丁继胜, 等. 多波束反向散射强度数据处理研究[J]. 海洋学报(中文版), 2006, 28(2): 51-55. Tang Qiuhua, Zhou Xinghua, Ding Jisheng, et al. Study on processing of multibeam backscatter data[J]. Acta Ocean ologica Sinica, 2006, 28(2): 51-55. |
[9] |
Beaudoin J, Hughes C, Ameele E, et al. Geometric and radiometric correction of multibeam backscatter derived from Reson 8101 systems[C]//Center for Coastal and Ocean Mapping. Canadian Hydrographic Conference. Toronto: Canadian Hydrographic Association, 2002: 242-266.
|
[10] |
阳凡林, 暴景阳, 胡兴树. 水下地形测量[M]. 武汉: 武汉大学出版社, 2016: 195-201. Yang Fanglin, Bao Jingyang, Hu Xingshu. Oceanic Surveying and Mapping[M]. Wuhan: Wuhan University Press, 2016: 195-201. |
[11] |
金绍华, 翟京生, 刘雁春, 等. Simrad EM多波束反向散射强度数据精处理研究[J]. 测绘科学, 2010, 35(2): 106-108. Jin Shaohua, Zhai Jingsheng, Liu Yanchun, et al. Study on processing of Simrad EM multibeam backscatter data[J]. Science of Surveying and Mapping, 2010, 35(2): 106-108. |
[12] |
Jensen F, Kuperman W, Porter M, et al. Computational Ocean Acoustics[M]. New York: Springer-Verlag New York, 2011: 73-81.
|
[13] |
Fisher F, Simmons V. Sound absorption in sea water[J]. The Journal of the Acoustical Society of America, 1977, 62(3): 558-564. |
[14] |
Urick R. The absorption of sound in suspensions of irregular particles[J]. The Journal of the Acoustical Society of America, 1948, 20(3): 283-289. |
[15] |
Richards S. The effect of temperature, pressure, and salinity on sound attenuation in turbid seawater[J]. The Journal of the Acoustical Society of America, 1998, 103(1): 205-211. |
[16] |
McMaster R. A statistical analysis of mathematical measures for linear simplification[J]. The American Cartographer, 1986, 13(2): 103-116. |
[17] |
周丰年, 赵建虎, 周才扬. 多波束测深系统最优声速公式的确定[J]. 台湾海峡, 2001, 20(4): 411-419. Zhou Fengnian, Zhao Jianhu, Zhou Caiyang. Determination of classic experiential sound speed formulae in multibeam echo sounding system[J]. Journal of Oceanography in Taiwan Strait, 2001, 20(4): 411-419. |
[18] |
Douglas D, Peucker T. Algorithms for the reduction of the number of points required to represent a digitized line or its caricature[J]. Cartographica, 1973, 10(2): 112-122. |