从零掌握Bellhop：二维水声射线仿真与Matlab自动化实践

1. 项目概述：从“听声辨位”到精准建模

在海洋这个庞大而复杂的声学环境中，声音的传播路径从来都不是一条简单的直线。无论是潜艇的隐蔽航行、海底资源的勘探，还是海洋环境的长期监测，都离不开对声波如何在水中“行走”的精确预测。这听起来有点像武侠小说里的“听声辨位”，但背后的科学要严谨和复杂得多。今天要聊的，就是一个在水声学研究和工程实践中堪称“瑞士军刀”的工具——Bellhop。这个项目，就是围绕Bellhop在二维海水声速环境下的声线仿真展开的。

简单来说，Bellhop是一个用于计算水下声场（特别是声线轨迹）的经典模型。它不像一些复杂的全波模型那样试图求解整个波动方程，而是采用了“射线声学”的近似。你可以把它想象成在海洋中发射无数条“声线”，然后追踪每一条声线在复杂声速剖面下的弯曲路径。这种方法计算效率高，物理图像直观，特别适合分析中高频声波在深海或变化平缓环境中的传播，比如分析声呐的作用距离、定位误差，或者设计水声通信系统的链路。

为什么是“二维”？在实际操作中，我们经常假设海洋环境在水平方向上是均匀的，或者只关心声源-接收器连线所在的垂直剖面内的传播情况。这样就把一个三维问题简化成了二维，极大地减少了计算量，同时抓住了声传播在垂直方向受声速剖面影响这一核心矛盾。对于初学者或需要快速评估方案的研究者来说，从二维仿真入手是绝佳的选择。

这篇文章，就是为你准备的。无论你是刚接触水声仿真的学生，还是需要在项目中快速验证声场特性的工程师，我将带你从零开始，一步步搭建一个基于Bellhop的二维声线仿真环境。我会详细拆解每一个输入文件的意义，分享如何用Matlab来驱动Bellhop并可视化结果，更重要的是，我会分享那些官方手册里不会写的“踩坑”经验和参数调校技巧。我们的目标很明确：让你不仅能“跑通”仿真，更能“读懂”仿真结果背后的物理意义，并灵活应用于自己的实际问题中。

2. Bellhop模型核心原理与输入文件拆解

要玩转Bellhop，首先得理解它的“游戏规则”。Bellhop本质上是一个根据给定的环境参数（水深、海底性质、声速剖面）和计算参数，求解声线方程的程序。它的输入全靠几个精心编排的文本文件，输出则是声线轨迹、传播损失等数据。理解每个文件里每个参数的含义，是避免“垃圾进，垃圾出”的关键。

2.1 环境文件（.env）：定义声波传播的舞台

.env文件是Bellhop的“总剧本”，它定义了仿真所处的物理环境。一个典型的二维仿真.env文件结构如下：

'My Bellhop Simulation' ! 第一行：标题，可任意填写，用于标识 1250.0 ! 第二行：频率（Hz），这是最重要的参数之一，直接影响声波波长和传播特性 1 ! 第三行：介质类型选项。1代表海洋水体，是最常用的设置 'SVF' ! 第四行：声速剖面类型。'SVF'表示由后面的离散点定义声速-深度关系 51 0.0 5000.0 ! 第五行：声速剖面点数、最小深度(m)、最大深度(m)。这里表示从海面0米到海底5000米，用51个点来描述声速变化。 0.0 1530.0 ! 声速剖面数据行。每行：深度(m) 声速(m/s)。必须从最小深度排到最大深度。 10.0 1525.0 20.0 1520.0 ... ... ! 中间点数据，描述了声速随深度的变化，比如温跃层导致的声速极小层 5000.0 1550.0 'A' 0.0 ! 海底边界类型和参数。'A'表示半无限液体海底，后面跟海底声速(m/s)。这是最简单的海底模型。 1600.0 0.0 1.0 0.0 ! 海底声速、密度(g/cm³)、衰减(dB/λ)、衰减幂指数。注意密度单位是g/cm³，1600.0就是1.6 g/cm³。 'R' ! 海面边界类型。'R'表示绝对刚性海面，声波全反射。这是最常用的简化。 0.0 ! 海面反射损失参数（对于'R'型，此参数无效，但必须保留一行）。 5000.0 ! 水深(m)。必须与声速剖面最大深度一致或小于它。 1 ! 声源数 0.0 100.0 ! 声源#1的深度(m)和指向性类型（100.0表示全向点源）。 1 ! 接收器数 0.0 5000.0 5001 ! 接收器#1的起始深度(m)、结束深度(m)、深度点数。这里表示在0-5000米水深均匀布置5001个接收器。

关键参数解析与避坑指南：

1. 频率（Freq）：这是仿真的“尺子”。频率越高，波长越短，声波更容易被小尺度不均匀性散射，射线理论的近似效果可能变差。通常，Bellhop适用于频率不太高（例如，对于海洋中尺度结构，几百Hz到几十kHz比较可靠）、环境变化不太剧烈的情况。如果你的频率设得过高（比如MHz级），得到的声线图可能会异常混乱，这时需要考虑更复杂的模型。
2. 声速剖面（SVF）：这是影响声线弯曲的“导演”。声速随深度变化（通常由温度、盐度、压力决定）是导致声线弯曲的根本原因。输入的点数要足够多，以平滑地描述声速变化，特别是在声速梯度变化剧烈的层（如温跃层）。点数太少会导致声线轨迹计算不准确。一个经验法则是，在声速变化剧烈的深度区间，采样间隔要小（如1-5米），在变化平缓的深海，间隔可以大一些（如50-100米）。
3. 海底边界：这是最容易出错的“配角”。'A'（半无限液体）模型最简单，但仅适用于声波几乎不穿透的硬质海底。对于松软的沉积层，需要使用'F'（流体海底）或'S'（弹性固体海底）模型，并需要额外提供海底层的厚度、声速、密度、衰减等参数。如果海底模型设置不当，会导致计算出的海底反射损失和折射行为完全错误，从而严重影响远距离传播损失的预测。
4. 接收器设置：接收器的深度范围和点数决定了输出声场图的“分辨率”。深度点数越多，图越精细，但计算量也越大。需要根据关注的重点深度区域（如声源深度附近、海底反射区）来合理设置。对于初步分析，可以设置较稀疏的点；对于精细分析，则需要加密。

注意：.env文件对格式极其敏感！多余的空格、错位的行、注释符!的位置不对，都可能导致Bellhop读取失败。最稳妥的方法是先从一个能成功运行的例子文件开始修改。

2.2 声线文件（.ray）：设定追踪的“探针”

.ray文件告诉Bellhop要发射哪些声线，以及如何发射。它控制着仿真的“视角”。

'Rays for My Simulation' ! 标题 51 ! 声线数目 -20.0 20.0 ! 声线发射角度的下限和上限（度，相对于水平面）。-20°表示向下发射，+20°表示向上发射。 0.0 ! 声源深度（m），通常应与.env文件中声源深度一致。 0.0 5000.0 5001 ! 声线追踪的深度范围（m），通常覆盖整个水深。 10000.0 ! 声线追踪的最大距离（m）。声线超出此距离或触及边界后停止追踪。

核心选择与经验：

• 声线数目与角度范围：声线数目决定了声场“采样”的密度。数目太少，声线图看起来稀疏，可能漏掉重要的传播路径；数目太多，计算时间变长。一个实用的方法是先使用中等数量（如51条）进行试算，观察声线覆盖情况，再进行调整。角度范围需要覆盖所有可能到达接收区域的声线。例如，如果声源在浅海，主要能量可能集中在靠近水平的方向；如果声源在深海，向上和向下发射的声线都可能很重要。
• 最大距离：这个值需要根据你关心的传播距离来设定。设得太小，可能看不到完整的会聚区或反转点；设得太大，又会浪费计算资源。通常可以先设一个较大的值，根据第一次仿真结果再调整。

2.3 其他输入文件：控制计算与输出

除了上述两个核心文件，通常还会需要：

• .bty文件：定义海底地形起伏。在二维仿真中，如果海底不是平的，就需要这个文件。它给出了距离-深度对，描述海底的轮廓。
• .trc文件：定义海面起伏（波浪）。用于考虑动态海面对声传播的影响，属于更高级的用法。
• .ati文件：定义海底声学属性横向变化。用于非均匀海底环境。

对于入门和大多数应用，.env和.ray文件就足够了。.bty文件在涉及斜坡、海山等地形时非常有用。

3. 基于Matlab的自动化仿真流程搭建

手动编写和修改这些文本文件非常繁琐，且容易出错。利用Matlab进行自动化驱动是提升效率的最佳实践。下面我将构建一个完整的Matlab脚本框架，实现从参数设置、文件生成、调用Bellhop到结果可视化的全流程。

3.1 环境参数设置与文件生成函数

首先，我们创建一个函数来生成.env文件。这样做的好处是参数集中管理，修改方便。

function writeBellhopEnv(filename, titleStr, freq, depth, c_profile, bottom_opt, src_depth, rcvr_depth_range, rcvr_num) % WRITEBELLHOPENV 生成Bellhop环境文件(.env) % 输入参数： % filename - 输出文件名（不含后缀） % titleStr - 仿真标题 % freq - 频率 (Hz) % depth - 水深 (m) % c_profile - Nx2矩阵，[深度(m), 声速(m/s)] % bottom_opt - 结构体，包含海底参数 % .type - 海底类型，如 'A' % .c - 海底声速 (m/s) % .rho - 海底密度 (g/cm^3) % .att - 衰减 (dB/λ) % src_depth - 声源深度 (m) % rcvr_depth_range - 接收器深度范围 [z_min, z_max] (m) % rcvr_num - 接收器数量 fid = fopen([filename, '.env'], 'w'); % 写入标题和频率 fprintf(fid, '''%s'' ! TITLE\n', titleStr); fprintf(fid, '%.1f ! FREQ (Hz)\n', freq); fprintf(fid, '1 ! NMEDIA\n'); fprintf(fid, '''SVF'' ! SSPOPTION\n'); % 写入声速剖面 [npts, ~] = size(c_profile); fprintf(fid, '%d %.1f %.1f ! NPTS, ZMIN, ZMAX\n', npts, min(c_profile(:,1)), max(c_profile(:,1))); for i = 1:npts fprintf(fid, '%.1f %.1f\n', c_profile(i,1), c_profile(i,2)); end % 写入海底参数 fprintf(fid, '''%s'' %.1f ! BOTTOM TYPE & C\n', bottom_opt.type, bottom_opt.c); fprintf(fid, '%.1f %.2f %.2f 0.0 ! BOTTOM: C, RHO, ATTEN, ATTEN_POWER\n', ... bottom_opt.c, bottom_opt.rho, bottom_opt.att); % 写入海面参数（假设为刚性） fprintf(fid, '''R'' ! SURFACE TYPE\n'); fprintf(fid, '0.0 ! SURFACE PARAM\n'); % 写入水深 fprintf(fid, '%.1f ! DEPTH (m)\n', depth); % 写入声源 fprintf(fid, '1 ! NSD\n'); fprintf(fid, '%.1f 100.0 ! SD, SRC_TYPE (100=omni)\n', src_depth); % 写入接收器 fprintf(fid, '1 ! NRD\n'); fprintf(fid, '%.1f %.1f %d ! RD_RANGE, NR\n', rcvr_depth_range(1), rcvr_depth_range(2), rcvr_num); fclose(fid); disp(['环境文件已生成: ', filename, '.env']); end

同理，我们可以创建生成.ray文件的函数writeBellhopRay。接下来，在主脚本中组织整个仿真：

%% Bellhop 2D 声线仿真主脚本 clear; close all; clc; % 步骤1：定义仿真参数 simTitle = 'DeepWater_Example'; freq = 1000; % 1 kHz waterDepth = 5000; % 5 km % 构建一个典型的深海声速剖面（存在声速极小层） z = (0:10:waterDepth)'; c = 1500 + 0.1*z - 0.001*(z-1000).^2; % 一个简单的抛物线模型，在1000米附近形成极小值 c_profile = [z, c]; % 设置海底参数（半无限液体） bottom.type = 'A'; bottom.c = 1600; % m/s bottom.rho = 1.8; % g/cm^3 bottom.att = 0.5; % dB/λ % 声源与接收器 srcDepth = 500; % 声源位于500米深度 rcvrDepthRange = [0, waterDepth]; rcvrNum = 501; % 接收器数量，决定传播损失图垂直分辨率 % 声线参数 numRays = 101; rayAngleLimits = [-30, 30]; % 度 maxRange = 100e3; % 最大追踪距离 100 km % 步骤2：生成输入文件 envFile = simTitle; writeBellhopEnv(envFile, simTitle, freq, waterDepth, c_profile, bottom, srcDepth, rcvrDepthRange, rcvrNum); writeBellhopRay([envFile, '.ray'], simTitle, numRays, rayAngleLimits, srcDepth, [0, waterDepth], maxRange); % 步骤3：调用Bellhop可执行文件 % 假设bellhop.exe位于当前目录或系统路径 bellhopExe = 'bellhop.exe'; % Windows系统 systemCmd = sprintf('%s %s', bellhopExe, envFile); disp(['执行命令: ', systemCmd]); [status, cmdout] = system(systemCmd); if status ~= 0 warning('Bellhop执行可能出错，请检查输入文件。'); disp(cmdout); end % 步骤4：读取并可视化结果 % 读取声线轨迹文件 (.ray) rayData = readBellhopRay([envFile, '.ray']); % 读取传播损失文件 (.shd) - 这里需要专门的读取函数，例如从AcTUP工具箱或自编 % [PlotTitle, PlotType, freq, atten, Pos, pressure ] = read_shd( [envFile, '.shd'] ); % 由于.read_shd函数通常需要特定工具箱，此处先聚焦声线可视化。

这里需要自己实现或获取readBellhopRay函数来解析Bellhop输出的.ray文件。该文件格式是固定的，可以按照官方文档编写读取代码。

3.2 声线轨迹可视化与解读

成功运行Bellhop后，我们会得到.ray文件，里面存储了每条声线的坐标（距离，深度）。可视化是理解结果的第一步，也是最关键的一步。

function plotRayTraces(rayData, waterDepth, maxRange, srcDepth) % PLOTRAYTRACES 绘制声线轨迹图 % rayData: 结构体，包含每条声线的r（距离）和z（深度）数组 % waterDepth: 水深 % maxRange: 绘图距离范围 % srcDepth: 声源深度 figure('Position', [100, 100, 800, 400]); hold on; grid on; box on; % 绘制每条声线 for i = 1:length(rayData) r = rayData(i).r / 1000; % 转换为km z = rayData(i).z; plot(r, z, 'b-', 'LineWidth', 0.5); end % 标记声源 plot(0, srcDepth, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r'); % 绘制海底和海面 plot([0, maxRange/1000], [waterDepth, waterDepth], 'k-', 'LineWidth', 2); % 海底 plot([0, maxRange/1000], [0, 0], 'k-', 'LineWidth', 2); % 海面 xlabel('距离 (km)'); ylabel('深度 (m)'); set(gca, 'YDir', 'reverse'); % 深度向下增加，这是水声学惯例 ylim([0, waterDepth*1.05]); xlim([0, maxRange/1000]); title('Bellhop 声线轨迹仿真'); legend('声线', '声源', 'Location', 'best'); hold off; end

从声线图中能看出什么？

1. 声速剖面的影响：如果声速随深度减小（表面声道），声线会向上弯曲；如果存在声速极小层（深海声道），声线会被“捕获”在该层上下摆动，传播距离极远。从图中声线的弯曲方向可以直接看出声速梯度的正负。
2. 影区（Shadow Zone）：由于声线弯曲，某些区域可能没有直达声线到达，形成“声影区”。在图中表现为某个距离和深度范围内没有声线穿过。这对于潜艇隐蔽和声呐探测至关重要。
3. 会聚区（Convergence Zone）：在深海，由于声线的周期性聚焦，会在特定距离上形成高强度的区域，即会聚区。在声线图上表现为多条声线在某个距离点交汇。
4. 海底/海面反射：观察声线与边界（海面、海底）的交互。刚性边界反射角等于入射角。如果设置了有损耗的海底，每次反射后声线强度会衰减。

实操心得：第一次看到声线图可能会觉得杂乱。一个很好的分析方法是，分角度区间观察。例如，单独绘制-5°到5°的声线，看看近水平发射的声线如何传播；再绘制-20°到-10°的声线，看看向下发射的声线如何与海底作用。这能帮你理清不同传播路径的贡献。

4. 传播损失计算与声场图绘制

声线轨迹给了我们路径信息，但最终我们往往更关心的是声能量随距离和深度的分布，即传播损失（Transmission Loss, TL）。Bellhop可以通过积分声线携带的强度来计算TL，并输出到.shd文件。

4.1 理解传播损失输出

.shd文件是一个二进制文件，存储了复声压场。传播损失TL通常定义为：TL = -20 * log10( |p| / |p0| )，其中p是接收点处的声压，p0是距离声源1米处的参考声压。

读取.shd文件需要专门的工具。Matlab社区有一个广为人知的工具箱叫AcTUP（Acoustic Toolbox User-interface & Post-processor），它包含了read_shd等函数，可以方便地读取Bellhop、KRAKEN等模型的输出。你可以从相关学术网站获取。

假设我们已经用read_shd函数读出了数据：

% 假设已使用AcTUP工具箱的read_shd函数 [PlotTitle, PlotType, freq, atten, Pos, pressure] = read_shd( [envFile, '.shd'] ); % Pos结构体包含了距离和深度轴信息 r_km = Pos.r.r / 1000; % 距离轴，单位km z = Pos.r.z; % 深度轴，单位m % pressure是复声压矩阵，尺寸为 (深度点数 × 距离点数) TL = -20 * log10( abs(pressure) ); % 计算传播损失 % 注意：这里计算的是相对损失。通常需要减去一个参考值，使得在声源处TL=0。 TL = TL - min(TL(:)); % 一个简单的归一化，使最小值为0

4.2 绘制二维传播损失等高线图

这是最直观的声场可视化方式。

figure('Position', [100, 100, 900, 400]); % 使用pcolor或contourf绘制 [R, Z] = meshgrid(r_km, z); contourf(R, Z, TL, 50, 'LineStyle', 'none'); % 50个颜色等级 colorbar; colormap(jet); % 或使用反转的灰度图 colormap(flipud(gray)) caxis([0 80]); % 设置颜色轴范围，0-80 dB是常见范围 xlabel('距离 (km)'); ylabel('深度 (m)'); set(gca, 'YDir', 'reverse'); title(['传播损失 (dB) - ', PlotTitle]); hold on; plot(0, srcDepth, 'w^', 'MarkerSize', 12, 'MarkerFaceColor', 'r'); % 标记声源 hold off;

解读声场图：

• 颜色深浅：代表传播损失大小。颜色越深（如蓝色），损失越大，声音越弱；颜色越浅（如红色/黄色），损失越小，声音越强。
• 明暗条纹：清晰的明暗相间条纹通常对应会聚区和影区。亮条纹表示声线聚焦，信号强；暗条纹表示声影区，信号弱。
• 垂直结构：靠近海面和海底，由于多次反射，通常会形成复杂的干涉条纹。
• 水平衰减：随着距离增加，整体颜色会变深，表示传播损失总体增大。

4.3 绘制特定深度的传播损失曲线

有时我们只关心某个特定深度（例如，水下航行器的工作深度）上，信号随距离如何衰减。

targetDepth = 1000; % 关注1000米深度 % 找到最接近目标深度的接收器索引 [~, idx] = min(abs(z - targetDepth)); tl_at_depth = TL(idx, :); figure; plot(r_km, tl_at_depth, 'b-', 'LineWidth', 2); grid on; xlabel('距离 (km)'); ylabel('传播损失 (dB)'); title(['在 ', num2str(targetDepth), ' 米深度处的传播损失']); % 标记会聚区位置（局部极小值点） [~, locs] = findpeaks(-tl_at_depth, 'MinPeakProminence', 5); % 找TL的谷值 hold on; plot(r_km(locs), tl_at_depth(locs), 'rv', 'MarkerFaceColor', 'r'); legend('传播损失', '会聚区', 'Location', 'best'); hold off;

这条曲线能直接告诉你，在目标深度上，每隔多远会出现一次信号增强（会聚区），以及信号随距离衰减的整体趋势，对于声呐作用距离评估和系统设计有直接指导意义。

5. 常见问题、调试技巧与模型局限性

即使按照步骤操作，你也可能会遇到各种问题。下面是我在无数次仿真中总结出的“避坑指南”。

5.1 仿真失败与错误排查

5.2 模型参数敏感性分析与校准

Bellhop的结果对某些输入参数非常敏感，了解这一点有助于合理解读结果。

1. 声速剖面：这是最敏感的参数。实测的声速剖面数据哪怕有小的误差，都可能导致预测的会聚区距离偏移几公里。务必使用尽可能准确的现场数据或可靠的海洋学模型输出。
2. 海底参数：对于远距离传播，海底反射损失是主要衰减源之一。海底声速bottom.c、密度bottom.rho和衰减bottom.att的猜测值会极大影响结果。如果有关注海底相互作用的项目，需要对海底参数进行反演或查阅区域地质调查资料。
3. 声线数目：声线数目不足会导致传播损失计算不准确，特别是在影区和会聚区边缘。一个简单的检验方法是：逐步增加声线数（如从51到101再到201），观察传播损失图是否收敛。当增加声线数后结果变化很小时，可以认为当前声线数已足够。

5.3 Bellhop模型的局限性认知

知道工具的边界和知道它的用法同样重要。

• 射线理论的局限：Bellhop基于高频近似（射线理论）。当波长与环境特征尺度（如内波、粗糙海面）相当时，衍射和散射效应变得重要，射线理论会失效。通常，频率低于1kHz或环境非常复杂时，需谨慎使用Bellhop，或将其结果与波动理论模型（如KRAKEN、SCOOTER）进行对比验证。
• 二维假设：忽略了水平方向的环境变化。如果实际环境中存在强烈的水平不均匀性（如锋面、涡旋），二维仿真结果可能与实际偏差较大。
• 确定性输入：Bellhop是确定性模型，输入确定的声速剖面和边界条件，得到确定的结果。它不包含海洋环境的随机起伏（如内波、湍流）带来的声场起伏。如需研究声场统计特性，需要结合随机介质理论或进行蒙特卡洛仿真。

一个实用的工作流建议：对于一个新的仿真场景，先从最简单的等声速（c_profile所有深度声速相同）和刚性海底开始。这时声线应该是直线，传播损失符合球面扩展规律（TL=20logR）。先确保这个基准案例能跑通且结果符合预期，然后再逐步引入真实的声速剖面和复杂海底模型。这种由简入繁的方法能帮你快速定位问题是出在环境设置上，还是出在工具使用本身。

最后，再分享一个调试时的小技巧：你可以让Bellhop输出更详细的信息。在.env文件末尾添加一行'A' *（注意星号前有空格），可以输出每个声线的振幅和相位信息，虽然文件会变大，但对于深度调试声线行为非常有帮助。仿真，尤其是水声仿真，是一个不断与物理现实和数值模型对话的过程。Bellhop提供了一个强大而直观的起点，但真正理解海洋如何“传递”声音，还需要你将仿真结果与实验数据、物理直觉以及其他模型反复对照。