简介这份资源面向水下声学、海洋探测与水下通信方向的研究者与工程人员提供 Bellhop 在海底地形起伏条件下模拟声波传播特性的 Matlab 源码适合具备一定射线声学基础、希望深入理解传播路径与衰减机制的中高级学习者。压缩包共 154 个文件约 33.63MB涵盖 m 源码脚本、prt 与 env 环境参数文件、bty 海底地形数据、ray 射线结果、shd 声场文件及 tif 图像等完整覆盖从参数配置到结果输出的计算链路。已有 1418 人学习下载具备较好的参考基础。读者可借助源码理解输入参数处理、射线追踪、散射计算与输出解析各环节结合声速剖面与地形数据开展定制化仿真并通过结果文件对照验证模型设置为水下通信信道评估、声纳性能预测等研究提供可复用的工具与排错思路。1. 从一张“声线图”说起Bellhop 到底能算什么第一次拿到这份 Bellhop 源码是在一个水声传播课题的联调现场。对方丢过来一个压缩包说“跑一下看看海底有山包的时候声场怎么变”。打开一看全是.m文件没有界面、没有按钮只有一堆函数和几个示例脚本。这正是 Bellhop 的典型形态它不是一个点开就能用的软件而是一套基于射线追踪Ray Tracing的水声传播计算工具用 Matlab 写成输入是环境文件输出是声线轨迹、传播损失、本征声线等结果。这份源码针对的场景很具体海底地形起伏条件下的传播特性。平坦海底的传播模型早就成熟了但真实海底有海山、有斜坡、有沟壑声线打上去会反射、会折射、会被遮挡传播损失曲线会出现明显的“阴影区”和“会聚区”偏移。Bellhop 用高斯波束Gaussian Beam方法处理这类问题比纯射线追踪更平滑能缓解焦散区能量发散的老毛病。适合谁做水声建模的研究生、做声呐性能预估的工程师、以及需要把传播模型嵌进更大仿真链路里的开发者。如果你只是想画一张漂亮的声线图交差它也能干但那就有点浪费了。2. 环境文件与运行链路一个算例从输入到出图2.1 目录里都有什么解压后先别急着run花两分钟把目录摸清楚后面改参数会省很多事。典型的 Bellhop Matlab 包结构大致是这样不同版本文件名略有出入但职责划分一致bellhop/ ├── bellhop.m % 主入口读环境文件、调度计算、写输出 ├── read_env.m % 解析 .env 环境文件 ├── ray2d.m / ray3d.m % 二维/三维射线追踪核心 ├── reflect.m % 海底/海面反射系数计算 ├── ssp.m % 声速剖面插值 ├── Arrival.m % 本征声线/到达结构提取 ├── examples/ % 算例环境文件.env 为主 │ ├── slope.env % 斜坡海底 │ ├── seamount.env % 海山地形 │ └── ... └── misc/ % 绘图、单位换算等辅助脚本.env文件是整个计算的“黑匣子入口”所有物理参数都写在这里。主脚本bellhop.m干的事很线性读环境 → 建声速场和边界 → 追踪射线 → 算场 → 落盘。理解这条链路比逐行读函数有用得多。2.2 环境文件长什么样Bellhop 的环境文件是纯文本分块书写块与块之间用单引号开头的注释行分隔。下面是一个带海底起伏的最小可用示例我按“能跑通”的标准精简过Seamount test ! 标题 50.0 ! 频率 Hz 1 ! 声源数 0.0 100.0 ! 声源深度范围 m 1 ! 接收器数 0.0 200.0 ! 接收深度范围 m 1 ! 接收距离数 0.0 5000.0 ! 接收距离范围 m CVW ! 边界条件C海面真空, V海底真空, W水柱 0.0 1500.0 0.0 ! 海面深度 声速 密度 100.0 1500.0 1.0 ! 水柱分层 200.0 1520.0 1.5 ! 海底界面起始 ALT ! 海底地形插值方式A解析, L线性, T样条 0.0 200.0 ! 地形点距离 深度 1000.0 180.0 2000.0 120.0 ! 海山隆起 3000.0 160.0 5000.0 200.0 A ! 声速剖面类型A解析 0.0 1500.0 0.0 ! 剖面参数几个参数值得单独说。第一行的频率决定波长波长又决定射线密度够不够频率越高需要的射线数越多否则会聚区算不准。CVW里的W表示水柱内是流体海底用V表示真空边界只是简化真实算例一般会换成A声学半空间并给出海底声速和衰减。地形块的ALT是插值方式T样条在起伏剧烈时容易过冲我一般先用L线性插值验证确认地形没读错再换样条。2.3 跑起来三条命令和一次检查Matlab 里进入包目录设置路径然后调用主函数% 添加路径注意 genpath 会把 examples 也加进来 addpath(genpath(pwd)); % 指定环境文件不带扩展名输出前缀同名 bellhop(examples/seamount); % 结果默认写到 seamount.shd传播损失和 seamount.ray声线bellhop的第一个参数是环境文件名不带.env后缀第二个可选参数是输出文件名前缀不传就沿用输入名。跑完先别急着画图用dir看一眼生成了哪些文件。常见输出有.shd复数声压场二进制、.ray射线坐标、.arr到达结构。如果只生成了.ray没有.shd多半是环境文件里接收器块写错了或者CVW那行的边界类型和后面的海底参数对不上。画传播损失最省事的办法是调包里的绘图脚本或者自己读.shd% 读取 shd 文件返回声压矩阵和距离/深度网格 [pressure, pos] read_shd(seamount.shd); tl -20*log10(abs(pressure) eps); % 传播损失 dB figure; pcolor(pos.r, pos.z, tl); shading interp; set(gca, YDir, reverse); % 深度向下 xlabel(距离 m); ylabel(深度 m); title(海底起伏条件下的传播损失);read_shd的返回维度通常是“深度 × 距离”pos.r和pos.z是对应的坐标向量。eps是为了防止对零取对数这个细节不加的话图上会出现-Inf白条第一次跑的人容易被吓到。3. 核心模块拆解射线追踪、反射与声场合成3.1 射线追踪主循环Bellhop 的射线追踪本质是解一组射线方程用龙格库塔类方法步进。核心逻辑在ray2d.m里简化后的骨架如下function ray trace_ray2d(env, theta0) % theta0: 初始掠射角单位弧度 % env: 已解析的环境结构体声速场、地形、边界 step env.step; % 步长默认与波长相关 r 0; z env.zs; % 起点距离 0深度为声源深度 phi theta0; % 当前射线方向角 ray [r, z]; % 轨迹缓存 while r env.rmax z env.zmin z env.zmax % 当前点的声速与声速梯度 c interp_ssp(env.ssp, z); dc gradient_ssp(env.ssp, z); % 射线方程dr/ds cos(phi)/c, dphi/ds -dc/dz * cos(phi)/c dr step * cos(phi) / c; dphi -step * dc * cos(phi) / c; r r dr; phi phi dphi; z z step * sin(phi); % 边界处理触底或触面则反射 if z interp_bottom(env.bathy, r) [phi, z] reflect_bottom(env, r, z, phi); elseif z 0 phi -phi; z 0; % 海面理想反射 end ray [ray; r, z]; end end这段代码有三个关键点。step不是越小越好太小会让长距离算例慢到不可接受太大则在地形陡变处直接穿模经验值是取波长的若干分之一同时在地形梯度大的区域自适应减半。dphi那一项里的dc是声速对深度的梯度正是它让声线在跃层附近弯折这也是声速剖面必须给准的原因。边界判断里interp_bottom对地形做插值如果地形点太稀疏射线会在两个采样点之间“跳”过海山表现为声线莫名其妙穿进海底——这是起伏地形算例最常见的翻车点。3.2 反射系数与海底建模射线打到海底不是简单镜面反射反射系数随掠射角变化还和海底声速、密度、衰减有关。reflect.m里通常用瑞利反射公式function R reflect_coeff(theta, c1, rho1, c2, rho2, alpha) % theta: 掠射角相对界面 % c1/rho1: 水介质声速/密度c2/rho2: 海底声速/密度 % alpha: 海底衰减 dB/波长 % 复声速把衰减折进去 k2 2*pi*f / (c2 * (1 1i*alpha/(40*pi*log10(exp(1))))); p sin(theta) * c1 / c2; % 斯涅尔定律 if p 1 R 1; % 全内反射 return; end cos2 sqrt(1 - p.^2); Z1 rho1 * c1 / sin(theta); Z2 rho2 * c2 / cos2; R (Z2 - Z1) / (Z2 Z1); endalpha的单位容易搞混有的文献用 dB/λ有的用 dB/m源码里一般按 dB/λ 处理换算时乘频率。p 1的分支对应全内反射此时反射系数模为 1相位仍会变如果直接返回实数 1 会丢掉相位信息导致干涉条纹算错。这个细节在平坦海底算例里不明显一到起伏地形、多次反射叠加时就暴露出来。3.3 高斯波束合成声场纯射线追踪在焦散区能量趋于无穷Bellhop 用高斯波束把每条射线的能量“抹开”再叠加得到平滑的声压场。合成逻辑大致是对每条射线沿轨迹计算波束宽度和曲率在接收点附近做高斯加权求和。function pressure beam_sum(rays, receivers, freq) % rays: 所有射线的轨迹与相位 % receivers: 接收点网格 [r, z] pressure zeros(size(receivers, 1), 1); k 2*pi*freq / 1500; % 参考波数 for i 1:numel(rays) ray rays(i); % 每条射线对每个接收点的贡献 for j 1:size(receivers, 1) dr receivers(j,1) - ray.r(end); dz receivers(j,2) - ray.z(end); % 高斯窗宽度由波束参数控制 w ray.beam_width; amp ray.amplitude * exp(-(dr^2 dz^2) / (2*w^2)); pressure(j) pressure(j) amp * exp(1i*k*ray.phase); end end endbeam_width是波束宽度它和频率、射线管面积有关源码里由射线间距和曲率推导。这个参数调小了声场退化成离散亮点调大了会聚区和阴影区的边界被抹平传播损失曲线失去细节。实际调参时我一般先固定射线数比如 2000 条扫一遍波束宽度看会聚区位置是否稳定稳定了再增加射线数提高精度。4. 避坑与排查起伏地形算例最容易栽的五个地方4.1 声线穿进海底现象声线图上射线直接穿过海底地形线继续在水下“飞行”。原因地形插值点太稀疏或者插值方式选了样条导致过冲射线步进时跨过了界面。解决把地形采样加密到与射线步长同量级插值方式先用L线性在ray2d.m的边界判断里加一步“半步回退”检测到越界就把步长减半重算。4.2 传播损失出现大片空白现象pcolor图上某些区域是白色或-Inf。原因该区域没有射线到达阴影区声压为零取对数后变-Inf也可能是接收网格没覆盖到。解决加eps防零是治标治本是增加射线数或放宽初始角度范围让波束能覆盖阴影区边缘。如果物理上确实是声影区空白就是正确结果别硬填。4.3 会聚区位置对不上文献现象算出来的会聚区距离比预期偏近或偏远。原因声速剖面单位或深度方向搞反了海底声速给成了淡水值频率和射线数不匹配。解决先画声速剖面确认深度向下递增检查.env里声速块的第一列是不是深度射线数按频率 × 距离 / 声速的量级估算不够就加。4.4 程序跑完没有任何输出文件现象命令窗口没报错目录里也没有.shd。原因环境文件名传错带了.env后缀输出路径没有写权限接收器块的距离范围写成了 0。解决bellhop的参数不带后缀用绝对路径试一次接收距离范围至少给一个非零值。4.5 三维算例内存爆掉现象跑ray3d时 Matlab 卡死或报内存不足。原因三维射线数按角度平方增长接收网格也是二维内存是二维算例的几十倍。解决先降射线数跑通流程确认环境文件无误再逐步加接收网格用稀疏采样出图时再插值必要时分距离段分批计算。5. 二次开发把 Bellhop 嵌进自己的仿真链路5.1 改造思路从“跑脚本”到“调函数”原始包是脚本式用法每次改参数都要动.env文件批量扫参很痛苦。二次开发第一步是把环境解析和计算拆开让参数从 Matlab 结构体传入而不是从文件读。做法是保留read_env.m作为文件入口另写一个env_from_struct.m把结构体转成内部格式主计算函数只认内部格式。这样批量跑 100 组海底坡度时不用生成 100 个.env文件。% 批量扫参示例海底坡度从 2° 到 10° slopes 2:2:10; results cell(numel(slopes), 1); for i 1:numel(slopes) env default_env(); % 默认参数结构体 env.bathy make_slope(slopes(i)); % 生成对应地形 env.freq 50; results{i} bellhop_core(env); % 纯计算不落盘 endbellhop_core是我建议新增的中间层输入结构体输出声压矩阵和射线不碰文件系统。落盘和绘图留给外层脚本。这样单元测试也好写给定一组已知解析解的平坦海底参数断言传播损失在容差内。5.2 验证方法拿平坦海底当“后悔药”改完核心代码最怕的是“看起来能跑其实算错了”。我的习惯是保留一个平坦海底的基准算例它的传播损失有近似解析解柱面扩展加海底反射修正每次改动后先跑基准误差超过 1 dB 就回滚。起伏地形没有解析解就用“收敛性验证”射线数翻倍看会聚区位置和传播损失曲线是否稳定稳定了才认为结果可信。验证项方法通过标准平坦海底对比解析近似传播损失误差 1 dB射线收敛射线数 1000→2000→4000会聚区距离变化 1%能量守恒全角度积分声强与源级偏差 0.5 dB地形连续性地形加密一倍声线轨迹不跳变5.3 一个具体技巧把波束宽度做成自适应固定波束宽度在起伏地形下很别扭平坦段够用海山附近要么太粗要么太细。我一般按局部射线管面积做自适应射线管被地形压缩时减小波束宽度避免能量“糊”到阴影区里。实现上在beam_sum里把w从标量改成随射线索引变化的向量由射线间距和曲率实时算。改完之后会聚区边界明显更锐利代价是计算量增加约两成换来的是地形遮挡效应算得准。从那以后我每次改 Bellhop 的核心函数都强制先跑一遍平坦海底基准再碰起伏算例这个习惯帮我挡掉过好几次“参数单位写错却看起来正常”的坑。这份源码不大但把环境文件、射线追踪、反射和波束合成这几块吃透再按上面的思路拆出可调用的中间层它就能从一个算例脚本变成你仿真链路里的一块积木。希望帮到你。本文还有配套的精品资源点击获取