简介tvapm.zip 是一套面向水声学研究的 MATLAB 模型程序包旨在帮助声学研究人员、海洋工程技术人员与水下通信工程师理解声波在海洋环境中的传播规律以及散射、吸收、折射等复杂物理过程。压缩包以水声传播模型与海洋环境参数模拟为核心包含 TVFilter、OceanTVIR、prepare_bellhop、UnderProSim 等核心脚本可覆盖声速剖面设置、海面波浪生成、海底地形输入到声线追踪、信道滤波与目标收发位置的完整仿真链条。全包共 23 个文件以 18 个 m 源文件为主用于构建各类数值模型4 个 dat 数据文件提供声速、海底地形、信号等仿真输入1 份 pdf 说明手册则给出使用指导和模型背景整体约 406KB轻量易部署。当前已有 319 人学习使用。借助这批代码使用者可快速搭建水声信道仿真流程结合 Bellhop 射线模型开展声线追踪、海底散射与信号滤波分析并通过对照示例数据与可视化结果深入理解不同海洋环境参数对声传播的影响。包内脚本按功能拆分便于按需调用或二次开发适合具备一定 MATLAB 基础、希望将理论模型转化为可运行代码的水声研究者与工程师。1. 从 tvapm.zip 认识一个老水声模型包的用法水声学模型很少以“安装后双击运行”的形式存在更多是像 tvapm.zip 这样一个压缩包把 Fortran 源码、数据文件和输出样例全塞在一起在课题组之间互相转发。TVAPM 这个名字指向一个具体的声传播模型程序pocketi1i 则是压缩包里一个难以直读的子程序文件名通常和 Pocklington 积分方程相关负责弹性体散射或边界积分计算。它没有图形界面没有完整文档价值在于你能看到每一行数学近似也能用最少操作得到一条可信的传播损失曲线。这篇文章按拆包、编译、输入、验证的顺序讲清楚这类老水声模型包从代码到结果的全过程适合手头拿到压缩包却不知道从哪儿开始的人。2. 拆包识代码tvapm.zip 的文件构成与 TVAPM 模型类型判断2.1 先列清单再解压从文件名判断水声模型类型判断一个水声模型包属于哪一类不用急着读代码文件清单已经泄露了大半信息。射线模型、简正波模型、抛物方程模型和目标散射模型它们的源码命名习惯完全不同。tvapm.zip 里有ssp.dat、bottom.dat这类数据文件很常见关键看主程序对应的算法缩写出现在哪个文件名里。mkdir tvapm_src cd tvapm_src unzip ../tvapm.zip find . -type f | sortmkdir一个独立目录是为了避免几十个.f文件散落到当前目录find | sort排序后第一眼就能分清主程序、子程序和输入数据文件。如果压缩包内有子目录说明原开发者按模块组织过代码保留这个结构再编译会省很多事。模型类型典型算法/程序文件名线索射线模型Bellhop、Raysray、bounce、arrive、tloss简正波模型Kraken、Orcamode、field、kraken抛物方程模型RAM、Pade、TDPEpe、ram、pade、tdpe积分方程/散射模型Pocklington、BEMpocket、pock、scat、bem射线适合高频近程简正波适合低频浅海抛物方程适合距离相关环境积分方程则多用于目标散射。如果清单里出现多个pe开头的文件而几乎看不到模式相关文件优先按抛物方程模型处理。TVAPM 缩写里的 V 在不同代码里分别代表 velocity、variable 甚至 time-varying在没有读到主程序注释之前不要把这个缩写写进论文先确认代码里那句C ... MODEL到底怎么定义。2.2 pocketi1i 到底是什么不认识的子程序先看文件头pocketi1i 这类名字多半是压缩包处理造成的表面现象而不是源码本来的名字。它更像是 pocketi1 或 pocket_integral_1 被长文件名规则截断后的结果i1 表示第 1 个积分核对应 Pocklington 积分方程中的一维积分算子。要确认它的真实身份不需要解压全部文件直接看文件内容更快。unzip -p tvapm.zip pocketi1i | head -40unzip -p把文件内容写到标准输出而不落盘适合快速预览。如果第一行是C POCKLINGTON INTEGRAL EQUATION KERNEL那就直接确认了。如果输出是 Fortran 固定格式源码把它重定向成.f文件再参与编译即可。unzip -p没有任何输出时先回 2.1 的清单确认文件名大小写压缩包里PocketI1和pocketi1i在 Linux 下是两个字面不同的文件。2.3 用 READ 语句反推输入文件格式老水声模型没有 README 时输入文件格式的权威来源是主程序里的 READ 语句。OPEN语句列出每个文件单元号READ语句列出每个变量的排列顺序。grep -n OPEN tvapm.f | head -30 grep -n -A4 READ(10 tvapm.f第一条命令告诉你要准备哪些文件第二条命令告诉你单元号 10 对应的输入文件里的列结构。常见模式是READ(10,*) nz后紧跟READ(10,*) (z(i), c(i), i1,nz)这个文件就是声速剖面第一行是采样点数之后每行是深度和声速。*表格式读入对空格数量和空行不敏感但对列顺序敏感所以不要随意调整输入文件的列位置。提示改输入文件前先复制一份原始文件。老模型的读入格式经常和某个数组维度绑定改了行数忘记改第一行点数程序通常会读出一串巨大数然后崩溃。3. gfortran 编译 TVAPM水声模型 Fortran 代码的适配要点3.1 老代码编译前先看 4 个编译选项水声模型的 Fortran 代码大多写于上世纪九十年代前后按 FORTRAN 77 固定格式组织。gfortran 能编译它们但默认选项会拦住一部分“当年合法、现在太老”的写法。编译选项作用什么时候用-ffixed-line-length-none解除固定格式 72 列行宽限制遇到Line truncated报错时-fno-range-check关闭常量范围检查老代码里常数溢出导致编译失败时-finit-realnan未初始化实数变量自动设为 NaN运行结果每次都不一样时-fbacktrace -ffpe-trapinvalid,zero,overflow捕获浮点异常并打印行号运行到一半崩溃时先用file tvapm.f确认是ASCII Fortran program text再用head -20 tvapm.f看第一列是否出现C注释或行号。只要第一列有C它就是固定格式-ffixed-line-length-none基本是必加的。3.2 最小编译命令和自带算例没有 Makefile 时最省事的方式是让编译器一次性处理目录下全部源码。前提是先确认没有多余的测试文件混在里面ls *.f扫一眼再动手。cd tvapm_src gfortran -O2 -g -ffixed-line-length-none -o tvapm *.f -lm-O2是老代码常用的优化级别-g保留调试信息-lm链接数学库。*.f会被 shell 展开成全部固定格式源文件如果包内同时存在.f90自由格式文件需要把它们单独列出并追加到命令末尾不能混在*.f里。编译输出如果只有 warning可以继续跑如果有 error优先看第一个报错位置老代码的语法错误常常是连续 5 行里每行带一个续行符而固定格式根本不认。编译通过后不要立刻跑自己的环境先找包内自带的输入样例。输入文件通常是*.dat、*.inp或*.env运行一次自带算例确认输出文件的行列数符合预期再替换成自己的声速剖面。3.3 三个高频问题LAPACK、重复符号、浮点异常第一个高频错误是链接阶段报undefined reference to dgesv_或zgesv_说明模型依赖 LAPACK 线性代数库。安装开发库后把链接库追加在源文件后面重新编译sudo apt install liblapack-dev libblas-dev gfortran -O2 -g -ffixed-line-length-none -o tvapm *.f -llapack -lblas -lm第二个问题是老代码自带 LINPACK 子程序系统库又提供同名函数报multiple definition。优先尝试删掉源码里自带的求解子程序因为它们通常可以从 LAPACK 等价替换如果不敢删把源码里的重复子程序改名后再编译。第三个问题是运行期浮点异常通常是数组越界、除零或开方负数。用调试选项重新编译一次gfortran -O0 -g -fbacktrace -ffpe-trapinvalid,zero,overflow -o tvapm_dbg *.f -lm ./tvapm_dbg run.prm-O0关闭优化以保证行号对应关系-ffpe-trap会在第一个浮点错误处停止并打印调用栈。栈顶指向的源码行就是需要检查的地方。4. 水声学模型的输入参数声速剖面、海底边界与网格步长4.1 从 CTD 数据生成声速剖面输入文件声速剖面是水声模型里最敏感的输入直接决定折射路径和会聚区位置。CTD 仪器测到的是温度、盐度和深度要先换算成声速再写入模型输入文件。import numpy as np d np.loadtxt(ctd_raw.txt) depth d[:, 0] temp d[:, 1] sal d[:, 2] c (1449.2 4.6 * temp - 0.055 * temp**2 0.00029 * temp**3 (1.34 - 0.01 * temp) * (sal - 35.0) 0.016 * depth) np.savetxt(ssp.dat, np.column_stack([depth, c]), fmt%10.2f %10.4f)这段代码把 CTD 原始三列数据转成两列的ssp.dat第一列深度第二列声速。np.savetxt的fmt保证每个数至少占 10 个字符避免老程序用宽松格式读入时因为字段贴在一起出错。深度方向要检查单调性很多模型要求深度从 0 向下递增CTD 下放和回收过程混在一起时必须先排序去重。老程序在内部通常做线性插值或三次样条输入点太密反而容易让插值曲线产生毛刺。对 100 Hz 左右的低频传播声速剖面保留每 5 m 一个采样点就足够了深水区甚至可以用每 50 m 一个点。不要直接用未经平滑的 CTD 微结构数据去跑远距离传播那会在输出里制造大量假干涉条纹。4.2 海底边界参数怎么给定量化海底对低频声传播的影响远大于多数人的直觉。老模型通常把海底简化成流体半空间只要三个参数纵波声速cp、密度rho和衰减alpha。海底类型纵波声速 cp (m/s)密度 rho (g/cm3)衰减 alpha (dB/波长)软泥1450–16001.1–1.50.2–0.5细沙1650–18001.8–2.00.5–1.0粗砂/砾石1800–21002.0–2.20.8–1.5岩石大于 2500大于 2.40.1–0.3设置参数前先在源码里确认衰减单位。有的程序用 dB/波长有的用 dB/m两者相差频率相关倍数。用 grep 找注释是最快的办法grep -in alpha\|atten\|db tvapm.f | head -20如果注释里写ALPHA IN DB PER WAVELENGTH而你的参考书给的是 dB/m换算公式是dB/m dB/波长 × f / c其中f是频率Hzc是海底声速m/s。同一组海底底质数据50 Hz 和 500 Hz 下的等效吸收可能差 10 倍换频率时不要只改频率参数海底衰减也要跟着重算。4.3 频率、网格步长与收敛性检查水声模型里频率决定了一切尺度。网格步长必须能分辨水中波长典型参考如下参数100 Hz 浅海1000 Hz 近程水中波长c1500 m/s15 m1.5 m垂直网格 dz 参考1.5–3 m0.15–0.3 m水平步长 dr 初始值5–10 m0.5–1.0 m垂直网格 dz 取波长的 1/10 到 1/5水平步长 dr 的保守初始值取波长的 1/4。老式抛物方程模型对 dr 的限制比新算法更严格不要直接拿论文里的 10 m 步长套到你的代码上除非你确认它用的是 split-step Padé 这类允许大步长的方法。判断网格是否足够最可靠的做法是步长减半对比。dr从 10 m 减到 5 mdz从 2 m 减到 1 m各跑一遍两个结果在相同距离点上的传播损失差小于 0.2 dB就认为网格收敛。如果差异超过 0.5 dB优先减小 dz垂直分辨率不足造成的相位误差远比为省时间跑粗网格严重。5. 用 TL 曲线和能量守恒验证 tvapm.zip 跑出的水声模型结果5.1 从二维输出中提取指定接收深度的 TL水声模型输出通常是一张二维表第一列是距离后面每一列对应一个接收深度。先head -20看表头确定目标接收深度在第几列然后用 awk 提取。awk NR1 $10 {print $1, $4} tl.out tl_rx.txt如果接收深度为 30 m垂直网格 dz 为 2 m对应的列号大约是30/2 2 17也就是第二个字段改成$17。NR1跳过表头$10过滤掉距离为 0 的边界行。提出来之后用 gnuplot 快速看一眼曲线形态plot tl_rx.txt using 1:2 with lines正常结果应该是光滑的衰减曲线叠加规则的低频干涉条纹。如果曲线出现锯齿状跳变先查声速剖面有没有深度不连续再查海底参数单位是否统一这个顺序可以排除掉大部分输入文件错误。5.2 两层解析对照近场和远场斜率均匀无边界介质中自由场传播损失是TL 20*log10(r) alpha*r其中r是距离alpha是吸收系数。把模型跑在均匀声速、无海底的测试配置里近场曲线如果明显偏离这个公式说明程序内部能量归一化或初值设置有问题。浅海环境看远场斜率。Pekeris 波导远场传播损失逐渐从球面扩展的 20log10(r) 过渡到接近 10log10(r)曲线斜率变缓代表能量被波导束缚。实测 TL 曲线远场比计算值陡很多时优先怀疑海底衰减参数设置过大而不是改声速剖面。awk NR1 {print $1, $2 20*log($1)/log(10)} tl_rx.txt | headawk 里log是自然对数log($1)/log(10)得到常用对数。这条命令把球面扩展项叠加到输出上检查余量是否接近常数。余量持续快速下降说明模型里有额外损耗来源。5.3 收敛性检查清单和一个海底敏感性技巧拿到一组可用结果后按三个顺序检查第一grep -i nan\|inf tl.out看输出里有没有非数值第二dr 减半后与原始结果的 TL 差是否小于 0.2 dB第三声源深度改变 0.5 m 后近场 1 km 处结果是否出现剧烈变化如果变化超过几分贝说明某个格林函数或镜像源实现里存在深度奇点。最后做一个海底敏感性粗检把海底吸收 alpha 分别设为 0 和 1 dB/波长跑两遍提取同一接收深度的 TL 曲线记录 10 km、50 km、100 km 三个距离上的差值。差值小于 3 dB说明该频段能量主要由声速剖面控制不需要花大力气标定海底差值大于 10 dB说明海底反演才是模型精度的主要矛盾。把这个差值表存成sensitivity_bottom.txt放到算例目录里下次调环境参数时先读它而不是重新猜。本文还有配套的精品资源点击获取
