企业官网建设流程全解析

简介：面波处理及剖面连接软件是一套面向地震勘探与地球物理研究人员的专业工具，集合了面波处理软件CCSWSWIN与速度分层软件CCSWSMAP，覆盖从地震记录去噪、相位解缠、频域与时频分析到速度分层建模、层析成像及多测线剖面连接等关键环节，可用于工程地质评价和资源勘探。压缩包共91个文件，仅1.8MB，以exe主程序、htm联机帮助、gif/bmp图示、doc使用说明以及dll/sys/vxd驱动和加密狗支持文件为主，辅以少量源码与工程文件，便于安装调试与对照学习。目前已有918人学习下载。借助这两个模块，使用者可以建立地下速度模型、拼接连续剖面并验证模型准确性；doc文档覆盖软件操作与仪器联接，gif图解和htm页面则直观展示处理流程，适合从入门到进阶的物探技术人员快速掌握。

1. 面波处理及剖面连接软件：把人工拾取变成自动追踪，省下的是实打实的半天工作量

在工程物探圈子里，面波处理一直带点“玄学”色彩——同一道野外数据，换两个人手工拾取频散曲线，结果能相差两成。这款面波处理及剖面连接软件做了一件反直觉的事：让频散曲线提取从人工鼠标点选变成算法自动峰值追踪，再让单点一维横波速度剖面按测线坐标自动拼成二维剖面，人只需要在开头定参数、在结尾复核结果，中间环节全部由程序代劳。

这套软件面向覆盖层调查、路基评估、空洞探测等工程物探场景，适合三类人：一线物探工程师想把野外面波数据批量变成可交付的剖面图；研究浅地表结构的学生想完整复现相移法、频散反演和空间插值的算法链路；做科学计算软件的开发者想参考“数据读取—谱分析—反演—剖面拼接”的分层实现。核心处理链路分五段，每一段的中间结果都有图像可以肉眼检查，不是个跑完才知道对错的黑匣子。

2. 系统架构与技术选型：SEG-Y读取、频散处理、剖面连接的分层管线

2.1 技术栈为什么选 Python

面波处理在商业软件里通常是 C++/Qt 的天下，因为交互响应快、原生界面性能好。但这套项目选择了 Python 3.8 + NumPy/SciPy + PyQt5，原因是常规测线的数据体量根本撑不起“必须上编译语言”的负担：24 道检波器、每道几千个采样点，一次 FFT 加相移叠加在 NumPy 里就是几十毫秒的事，用户点完“处理”到频散谱图刷新出来，几乎察觉不到延迟。性能不构成瓶颈之后，代码可读性和改造成本就成了第一优先级，Python 明显更合适。

界面层用 PyQt5，绘图部分用 Matplotlib 的 FigureCanvas 嵌进 Qt 窗体，频散谱、一维反演曲线、二维剖面共用同一套绘图逻辑，输出 PNG、SVG 也是顺手的事。算法层完全依赖 NumPy 和 SciPy，FFT 用 numpy.fft，滤波用 scipy.signal，插值用 numpy.interp，没有引入 ObsPy 这种大而全的地震学库。不用的原因很实际：ObsPy 安装链路长、依赖一堆，很多工区现场的 Windows 电脑装上就是灾难；而 SEG-Y 读取本身不复杂，项目里手写了一个约 80 行的解析器，反而更容易按需修改。

反演部分没有用 scipy.optimize 里的梯度类算法。频散曲线反演的目标函数是层状介质正演算出的理论相速度与实测值的误差，这个函数在参数空间里多峰，梯度法很容易停在局部极小，这在反演任务里是出了名的坑。常见做法是改用启发式全局搜索，项目里选了粒子群（PSO），因为实现简单、调参只有四个参数，还不用求导数。

配置格式固定用 JSON，是刻意保持简单的决定。JSON 在 Python 里就是原生 dict，界面上的参数表单可以直接映射成控件，省掉一个解析层。有人建议换 YAML，但代码里有四个模块都在用 json.load，真要换会动到不止一处，收益也不明显。对一个单人维护的项目来说，省心比炫技重要。

2.2 目录结构：四个层级的工程布局

拿到源码第一件事不是跑 main.py，而是先看目录骨架。项目按“界面层、算法层、配置层、测试层”四层组织，核心算法不依赖 Qt，界面只做交互与展示。

surface_wave_processor/ ├── main.py # 程序入口，支持 single / batch / gui 三种模式 ├── requirements.txt ├── config/ │ ├── default_config.json # 默认处理参数 │ └── demo_config.json # demo 数据配套参数 ├── core/ # 算法层，禁止 import PyQt5 │ ├── io/ │ │ ├── segy_reader.py # SEG-Y 格式解析 │ │ └── text_reader.py # 通用文本道集读取 │ ├── processing/ │ │ ├── preprocess.py # 滤波、增益、去均值 │ │ ├── dispersion.py # 相移法频散谱计算 │ │ ├── pick_curve.py # 频散曲线自动拾取 │ │ └── inversion.py # PSO 反演 │ ├── profiling/ │ │ ├── connect_profile.py # 剖面连接 │ │ └── layer_track.py # 层位自动追踪 │ └── utils/ │ └── logger.py # 日志与中间状态记录 ├── ui/ # 界面层，不写算法 │ ├── main_window.py │ ├── spectrum_widget.py # 频散谱绘制控件 │ └── profile_widget.py # 剖面绘制控件 ├── data/ │ ├── demo/ # 24 道 / 48 道 / 60 道三组数据 │ └── output/ └── tests/ ├── test_dispersion.py └── test_profile_connect.py

这个结构里有一条硬约束：core 目录下任何模块不允许 import PyQt5，ui 目录下任何模块不允许直接写算法运算。违反这条约束的重构一眼就能看出来。这样分层的好处是，核心模块可以在命令行下直接跑单元测试。我一般拿到代码会先执行pytest tests/ -v，确认算法层没问题再开界面，否则界面上点出来的报错既有 Qt 的堆栈又有 NumPy 的维度信息，排查起来是双重折磨。

core 内部又按职责拆了三段：io 只管格式解析，processing 管频散谱、拾取、反演，profiling 管一维剖面拼成二维剖面。utils 里放的日志模块看起来不起眼，实际作用很大。面波处理中间步骤多，几十个测点批量跑下来，哪个测点死在哪个环节，日志里一行就要能定位，不然一个个查文件会非常痛苦。

data/demo 目录里放了三组模拟数据，分别是 24 道、48 道和一组明显加了噪声的 60 道数据。demo_config.json 把道间距 2 米、采样间隔 0.0005 秒、偏移距 5 米这些采集参数都配好了，第一次试跑直接用它，省得自己去理解每个参数的含义。千万不要直接改 default_config.json，正确的做法是先复制出一份新文件再改，保留原始配置当参照，后面排查问题会省很多力气。

2.3 数据管线：一个测点从原始记录到一维剖面的流程

整个处理管线拆成六个阶段，每个阶段都有独立模块。单测点处理和整条测线批量处理复用同一套函数，区别只在外层循环的调度方式。各阶段的输入输出能直接用一张表说清楚。

批量处理时，main.py 按测线目录里的文件顺序遍历，但不会一口气跑完。每完成一个测点的预处理和频散谱计算，程序就往 output/intermediate_status.json 写入进度。假如第 8 个测点反演不收敛，你不用回头重跑前 7 个，用--resume参数让程序从第 8 个继续。实际工区经常有几十个测点，一次跑几十分钟是常态，没有断点续跑这个机制，中途崩一次能把人折磨到怀疑人生。

单测点的执行流程很有代表性，拿 demo 数据来说，跑完整条“读取—预处理—频散谱—拾取—反演”链路后，output 目录下会多出频散谱图、频散曲线 JSON 和反演剖面图三个文件。先确认这三样东西正常生成，确认环境没有问题，再进入批量模式。批量模式跑完后会额外生成一条二维剖面图，也就是软件名字里“剖面连接”真正出成果的环节。

2.4 界面层与算法层的协作方式

界面层没有直接调 dispersion.py 里的函数，而是利用 QThreadPool 和 QRunnable 把批量任务放到后台线程。用户在主窗口点“开始处理”，界面先收集参数打包成 dict，交给一个 Worker 对象，Worker 在子线程里调用 core 层的算法函数，每完成一个测点就发射一次进度信号，主线程收到信号后更新进度条。取消任务则用一个 threading.Event 标记，后台循环每处理完一个文件就检查一次标记，然后安全退出。

这样设计是血泪教训换来的：第一版直接在按钮的槽函数里写for shot in shots: process(shot)，一跑批量界面就白屏、系统弹“未响应”，macOS 上更是直接把窗口变灰色。后来把所有耗时计算全部移出主线程，界面才真正像桌面软件。QThreadPool 的线程数固定为 1，不让多道数据并行，因为面波处理里的 FFT 和矩阵运算内存占用不低，多测点并行会让内存翻倍，中间文件写入也会互相竞争。串行处理虽然慢一点，但稳定得多。

3. 核心模块拆解：相移法频散谱与剖面连接插值的代码实现

3.1 相移法频散谱：从波形到能量图的代码实现

频散谱是整个软件的承重墙。只有把原始记录变换到“频率-相速度”域，后续的自动拾取和反演才有输入。项目里采用的是相移法，它只需要做一次 FFT 加后续相移扫描，计算效率在各类频散谱算法里性价比最高。

import numpy as np def dispersion_spectrum(data, dt, dx, fmin, fmax, vmin, vmax, nv=200): """ 相移法计算频散谱 data: shape (n_receivers, n_samples) 的二维波形数组 dt : 采样间隔，单位秒 dx : 道间距，单位米 fmin, fmax : 关注频率范围，单位 Hz vmin, vmax : 相速度搜索范围，单位 m/s nv : 速度扫描点数 返回 freqs, velocities, disp_map disp_map shape 为 (nv, n_freq)，方便直接 imshow """ n_recv, n_samples = data.shape # 每一道做 FFT，得到频率域谱 spec = np.fft.fft(data, axis=1) freqs = np.fft.fftfreq(n_samples, d=dt) # 只保留有效正频率段 freq_mask = (freqs >= fmin) & (freqs <= fmax) freqs_sel = freqs[freq_mask] spec_sel = spec[:, freq_mask] velocities = np.linspace(vmin, vmax, nv) distances = np.arange(n_recv) * dx disp_map = np.zeros((len(freqs_sel), nv)) for i, f in enumerate(freqs_sel): # 相移矩阵：道数 x 速度数 phase = np.exp(-1j * 2.0 * np.pi * f * np.outer(distances, 1.0 / velocities)) # 对接收道做相移叠加 stacked = np.abs(spec_sel[:, i] @ phase) / n_recv # 每行做 0~1 归一化 stack_max = stacked.max() if stack_max > 1e-12: stacked /= stack_max disp_map[i, :] = stacked return freqs_sel, velocities, disp_map.T

逻辑说明：核心思想是给每道频谱乘一个与偏移距、测试速度有关的相位因子，然后在道方向求和。当测试速度恰好等于某个频率对应的真实相速度时，各道相位基本一致，叠加后能量最大；速度不对时，各道相位互相抵消，能量自然变低。外层循环遍历频率，内层用一次矩阵乘法扫完 200 个速度，比双重 for 循环快一个数量级。

参数说明：dt 和 dx 必须来自观测系统，采样率 0.5ms 就传 0.0005，写错一个数量级整个谱全糊。fmin 和 fmax 根据震源和检波器频带决定，锤击面波有效频带一般在 5~80Hz。fmin 设太低会把地脉动低频噪声也扫进来，fmax 设太高会触发空间假频。vmin 和 vmax 决定谱的纵轴覆盖范围，我一般先跑一个宽松范围比如 100~1000m/s，看能量团集中在哪儿再缩窄。

这个函数里最容易翻车的是np.fft.fftfreq(n_samples, d=dt)的 dt。dt 单位是秒，很多从 MATLAB 转过来的人习惯传毫秒，结果所有频率都放大了 1000 倍，频散谱上原本应该在 20Hz 的能量条带跑到了 20000Hz，怎么看怎么不对。这个问题在避坑章节会再展开。

3.2 频散曲线自动拾取与反演

频散谱算出来后，下一步是把它压缩成一条频散曲线。所谓拾取，就是为每个频率找能量峰值对应的相速度。这一步最大的干扰有两个：低频端噪声能量高，峰值不稳定；高频端基阶与高阶模式混叠，峰值可能在多个速度之间跳来跳去。

import numpy as np from scipy.signal import medfilt def pick_dispersion_curve(freqs, velocities, disp_map, threshold=0.6, kernel=5): """ 从归一化频散谱中自动拾取频散曲线 freqs: 频率数组; velocities: 相速度数组 disp_map: shape (n_freq, n_vel) threshold: 低于该能量视为噪声 kernel: 中值滤波窗口，必须为奇数 返回 picked_freq, picked_vel """ nfreq = disp_map.shape[0] picked_f = [] picked_v = [] for i in range(nfreq): row = disp_map[i, :] if row.max() < threshold: continue peak_idx = np.argmax(row) picked_f.append(freqs[i]) picked_v.append(velocities[peak_idx]) picked_v = np.array(picked_v) if len(picked_v) > kernel: # 中值滤波去掉单点跳变 picked_v = medfilt(picked_v, kernel_size=kernel) return np.array(picked_f), picked_v

逻辑说明：程序对每个频率单独找能量峰值，低于 threshold 的直接丢弃，避免把噪声当信号。拾取完再做一次中值滤波，把相邻频率点上的细小抖动抹平。threshold 取 0.6 是项目里的经验默认值，信噪比低的测点要降到 0.3 左右，否则低频端会丢很多点；但 threshold 太低，高频端又容易把高阶模式拾进来。我一般会在批量跑之前先打开频散谱图扫一遍，定一个整条测线统一的 threshold，避免每个测点单独调参导致剖面在测点连接处出现不连续。

反演部分用的是粒子群算法，目标函数是层状介质正演算出的理论相速度与实测频散曲线的均方根误差，粒子位置代表各层横波速度。项目里把正演封装在 inversion.py 的 forward 函数里，按层状模型算出理论频散，再把误差返回给 PSO 迭代。PSO 参数固定为：粒子数 80、迭代 200 轮、惯性权重 0.6、个体学习因子 1.4、群体学习因子 1.4。单个测点跑下来大概 10 秒左右，批量几十个测点几分钟可以跑完。

真正需要小心的是层厚度和层数的设定，这直接决定反演有没有解。浅层 30m 探测目标，层厚按 1m、2m、4m、8m、15m 递增设置，而不是等厚分层。面波对浅层最敏感，浅部必须有足够细的层去约束，深部地层对频散曲线响应弱，只能用越来越厚的层去代表。

3.3 剖面连接：线性插值与层位追踪的组合

一维反演完成后，每个测点得到一组“深度-横波速度”曲线。把这些曲线沿测线连成二维剖面，是软件名字里“剖面连接”的含义。项目里先做线性插值保证全局连续，再用层位追踪把插值结果的折线感磨平。

import numpy as np def connect_profiles(x_coords, vs_profiles, depth_profiles, n_cols=200, n_rows=100): """ x_coords: 测点沿测线 x 坐标，单位 m vs_profiles: 每测点反演得到的横向波速度数组列表 depth_profiles: 每测点对应的深度数组列表 n_cols, n_rows: 输出剖面的网格密度 返回 xs, zs, vs_grid """ x_min, x_max = x_coords.min(), x_coords.max() z_max = max(dp[-1] for dp in depth_profiles) xs = np.linspace(x_min, x_max, n_cols) zs = np.linspace(0, z_max, n_rows) vs_grid = np.zeros((n_rows, n_cols)) for j, x in enumerate(xs): # 找当前 x 左右两个测点索引 if x <= x_coords[0]: left = right = 0 elif x >= x_coords[-1]: left = right = len(x_coords) - 1 else: right = np.searchsorted(x_coords, x) left = right - 1 # 深度方向逐步插值 for i, z in enumerate(zs): v_left = np.interp(z, depth_profiles[left], vs_profiles[left]) v_right = np.interp(z, depth_profiles[right], vs_profiles[right]) # 按水平距离线性加权 dx_total = max(x_coords[right] - x_coords[left], 1e-6) t = (x - x_coords[left]) / dx_total vs_grid[i, j] = (1 - t) * v_left + t * v_right return xs, zs, vs_grid

逻辑说明：二维网格里每个点，先判断它落在哪两个测点之间，再在深度方向用 np.interp 分别取左右两个一维剖面的速度值，最后按水平距离加权。这个写法在测点间距大时会有明显折线，看起来像锯齿状层位。项目里后来又加了一个 layer_track.py 模块，先在前一步插值结果里搜索连续的速度等值线，再沿测线方向做中值滤波平滑，出图效果自然很多。

参数说明：n_cols 和 n_rows 是输出剖面纵向和横向网格密度，200x100 已经足够用于日常出图，不需要更密。测点间距超过 10m 时，线性插值会把中间可能存在的异常体直接抹掉，这时就要改用克吕金插值，把测点速度看作空间随机场，用变差函数建模后估计未知点。connect_profile.py 的接口已经预留了替换空间，第 6 章的二次开发会演示怎么把这段线性插值换成克里金。

4. 运行与配置：虚拟环境、JSON 参数和 Demo 数据验证

4.1 环境依赖与安装

项目依赖全部集中在 requirements.txt 里，包括 numpy、scipy、matplotlib、PyQt5、pytest。没有额外的大件依赖，这对工区电脑很友好。安装时推荐用 conda 建一个独立的 Python 3.8 虚拟环境，避免和系统 Python 混在一起。Python 3.9 和 3.10 也能跑，但 PyQt5 在 3.8 下验证最多，遇到奇怪的界面加载问题可以退回 3.8 试试。

conda create -n surfacewave python=3.8 conda activate surfacewave cd surface_wave_processor pip install -r requirements.txt

安装过程里最容易出问题的是 PyQt5 在离线环境下装不上。工区电脑如果没外网，提前把三个 whl 文件拷贝到现场，用pip install --no-index --find-links=./pkgs装离线包。另外 matplotlib 默认后端在无显示环境下会报错，如果通过远程终端跑批处理又不需要出图，可以在命令前加MPLBACKEND=Agg环境变量。

4.2 配置文件解读：一个 JSON 从采集参数控制到输出设置

项目所有处理参数都收敛到一个 JSON 文件里，好处是改动可追溯、方便复制。demo_config.json 打开后能看到几个明显的分段，每一段正好对应数据处理管线的一个环节。

{ "survey": { "line_name": "demo_line1", "receiver_count": 24, "receiver_spacing": 2.0, "shot_offset": 5.0 }, "preprocess": { "bandpass": [5, 12, 45, 80], "gain": "agc", "window_length_ms": 100 }, "dispersion": { "fmin": 5, "fmax": 80, "vmin": 100, "vmax": 1000, "nv": 200 }, "pick": { "threshold": 0.6, "medfilt_kernel": 5 }, "inversion": { "layer_thickness_m": [1, 2, 4, 8, 15], "pso_particles": 80, "pso_iterations": 200, "vs_range": [100, 800] }, "profile": { "x_step": 1.0, "z_step": 0.5, "interpolator": "linear" }, "output": { "output_dir": "data/output", "image_format": "png", "dpi": 100 } }

修改建议是：换一个工区的数据时，先复制一份 config 文件，把 line_name、receiver_spacing、fmin/fmax、vs_range 这几项按新工区采集参数改掉，其他保持默认。不要直接在默认配置上改，否则一次参数调乱了，下次想找回原始状态很难。

4.3 启动验证：单测点、批量和 GUI 三条路径

第一次跑项目，推荐先走单测点模式，用命令把处理链路完整跑一遍：

python main.py --mode single --config config/demo_config.json --shot data/demo/line1/shot001.dat

正常结束后 output 目录下应该出现三个文件：shot001_dispersion.png 是频散谱图，shot001_curve.json 是拾取出的频散曲线，shot001_inversion.png 是一维横波速度反演剖面。三样东西都在，说明环境没问题、算法模块能跑通。这时候再进批量模式：

python main.py --mode batch --config config/demo_config.json --survey data/demo/line1

批量模式会多生成一个 line1_profile.png，这就是最终连起来的二维剖面。如果批量跑到一半中断，带上--resume参数可以从断点继续，不需要重头再来。

GUI 模式验证更直观，命令是：

python main.py --mode gui

启动后从文件菜单打开 data/demo/line1/shot001.dat，界面左侧会显示原始波形，右侧是频散谱图。参数栏会自动加载 demo_config.json 里的值，点“开始处理”后流程和命令行完全一致，只是多了一个进度条。验证重点看频散谱上的能量条带：低频段能量团应该偏高速方向，比如 200~400m/s，高频段能量团往低速方向走，比如 100~200m/s，中间连成一条连续条带。如果看到的是均匀噪点，优先检查 dt 和 fmax 两个参数。

5. 避坑手册：面波处理软件开发的五个高频翻车点

5.1 频散谱能量图糊成一团

现象：demo 数据跑出来条带清晰，换成野外实测 24 道数据后，频散谱上没有明显能量集中，整幅图像均匀噪点。

原因：两个因素最可疑。一是空间假频，道间距 2m、道数 24 时，空间采样决定可恢复的最高频率有限，fmax 设到 80Hz 但高速层相速度在 600m/s 以上，能量发生混叠，谱图全乱；二是 dt 写错，采样间隔 0.0005 秒写成了 0.5 秒，FFT 频率轴整体放大 1000 倍。

解决：先用上限速度估算截止频率，常见做法是 f_max_limit = v_min / (2 * dx)，超过这个值的高频段直接截掉。然后检查 config 里 dt 有没有写对。这两个参数改好，80% 的“糊图”问题都能解决。

5.2 剖面连接出现不存在的台阶状断层

现象：测线上相邻两个测点之间，二维剖面在某一深度区间出现垂直速度突变，单看剖面像有断层，但地质资料显示该区域层位连续。

原因：两个测点反演时初始模型设置不一致，或拾取阈值不同导致频散曲线点数相差太大。面波反演是病态问题，同一组观测数据在不同初始模型下会收敛到不同局部极值，速度差可以达到 20% 以上，拼到同一张图上就是一个垂直台阶。

解决：全测线统一 threshold、medfilt_kernel、层厚度模型和 PSO 随机种子，把参数收敛到一个全局 config 里。检查时用 output/intermediate_status.json 看每个测点实际用的参数是否一致。如果只有一个测点异常，用--resume --shot 5单独重跑 5 号测点，对比频散曲线和反演拟合误差，必要时按区域经验收紧 vs_range。

5.3 SEG-Y 读取出错：头块字节位置对不上

现象：读某厂家的 SEG-Y 文件，程序直接报ValueError: unpack requires a buffer of 4 bytes，或者读出的道序完全乱掉。

原因：SEG-Y 格式有标准但实现五花八门，道头里采样点数和采样率的偏移量在不同采集软件里不一样，有的厂商还塞了私有字段。手写的解析器按标准偏移解析，遇到非标文件就会崩。

解决：在 segy_reader.py 里加一个“头块字段探测”逻辑，先读 240 字节道头，用几个候选偏移去解析采样点数，判断读出的值是否在合理范围（比如几百到几万），不在就试下一个偏移。每次解析把实际用到的偏移量打印到日志里，这样现场遇到新格式，看日志就知道该往哪个偏移修。

5.4 界面批量处理时假死

现象：点击“开始处理”后窗口立刻变“未响应”，Windows 上标题栏出现转圈，macOS 上直接灰屏。

原因：批量循环直接跑在了 Qt 主线程上。Qt 主线程要持续处理重绘和事件分发，CPU 密集型任务抢占线程后，窗口自然失去响应。这不是计算量大小的问题，哪怕只有十个测点，只要任务阻塞主线程几十秒，界面都会假死。

解决：把处理任务包装成 QRunnable，用 QThreadPool 提交到后台线程。做法是：

from PyQt5.QtCore import QThreadPool, QRunnable class ProcessWorker(QRunnable): def __init__(self, fn, *args): super().__init__() self.fn = fn self.args = args def run(self): # 后台线程里执行核心算法 self.fn(*self.args) pool = QThreadPool.globalInstance() pool.setMaxThreadCount(1) pool.start(ProcessWorker(process_survey, config_dict, event))

主线程只刷新进度条，取消按钮通过 threading.Event 通知 Worker 在文件边界退出。线程数设成 1，不要加多线程并行，面波处理里面 FFT 和矩阵运算本来就吃内存，多线程并行内存直接翻倍，中间文件写入还会互相干扰。

5.5 同一份数据跑两次反演结果差很多

现象：同一份数据同一天跑了两次反演，两条一维速度曲线相差 30% 以上，但两次日志里的处理参数完全一致。

原因：PSO 的粒子初始位置是随机生成的，早期版本没有固定随机种子，每次优化起点不同，结果自然不可复现。这在地球物理反演里尤其致命，因为反演本身病态，随机性会被放大。

解决：每次任务启动时用测点 ID 生成随机种子，例如rng = np.random.default_rng(int(shot_id.encode().hex(), 16) % 2**32)，并把随机种子、粒子数、迭代次数一起写进结果 JSON。从那以后我再跑批量任务前，都会强制检查 config 里有没有 random_seed 字段，没有就直接拒绝开始跑。

6. 二次开发指南：把剖面插值换成克里金并保持测试通过

6.1 替换插值内核：半小时接入克里金

connect_profile.py 里线性插值逻辑其实独立成了一个内部函数，替换起来手很快。常见的做法是在 core/profiling 下新增一个 kriging.py，实现空间变差函数建模和克吕金预测，然后在 connect_profile.py 里加一个开关字段，让 profile.interpolator 支持 "linear" 和 "kriging" 两个值。调用处不用大改，接口保持统一即可。

# core/profiling/connect_profile.py 内部伪代码 from core.profiling.kriging import kriging_interpolate if config["profile"]["interpolator"] == "kriging": vs_grid = kriging_interpolate( x_coords, vs_profiles, depth_profiles, variogram_model="spherical", max_neighbors=8 ) else: vs_grid = linear_interpolate( x_coords, vs_profiles, depth_profiles, n_cols=n_cols, n_rows=n_rows )

spherical 变差函数模型适合大多数浅表速度场，max_neighbors 限制每个估计点最多参与的测点数，设 8 表示只让邻近测点影响插值结果，避免远处测点把局部异常冲淡。改完后跑一遍 demo 数据，和原线性插值结果对比，正常情况下克里金结果更平滑，异常体边界更清楚。

6.2 新增一种数据格式导入器

工区偶尔会出现第三方采集仪器的自定义文本格式，这时候只需要在 core/io 下照 text_reader.py 的结构写一个 csv_reader.py，实现 parse_head 和 parse_trace 两个函数，再在 io/init.py 的格式注册表里加一行。读取器不需要做任何频散计算，职责边界保持干净，后续换格式不会影响算法层。

6.3 用回归测试确认没改坏

替换完插值算法后，跑pytest tests/test_profile_connect.py -k kriging，用 tests 里预置的合成层状模型验证：合成模型的真实速度场已知，插值结果与真实速度场的最大偏差应该小于设定阈值。最后再用 demo 数据出一张二维剖面图，和改动前的图并排检查边界连续性。我以前改完插值不跑测试直接拿工区数据出图，结果剖面拐出一个根本不存在的低速异常，交付后被审查单位打回重做。从那以后我给自己立了个规矩：任何算法改动，先跑单测，再用 demo 数据出图对照，最后才碰真实数据入库。希望帮到你。

本文还有配套的精品资源，点击获取