简介：这是一份面向地质探测与信号处理方向学习者的探地雷达数据处理程序源码包，适合具备一定编程基础、希望深入理解GPR数据解析流程的开发者与研究人员。程序围绕道内均衡、道间均衡、指数增益、频率域与时域滤波等核心环节展开，并配有二维、三维数据可视化功能，可用于优化雷达回波信号、提升深度剖面中地质特征的辨识度。压缩包共94个文件，约12.06MB，以h头文件与cpp源文件为主体，辅以obj编译中间文件、bmp与jpg图像、sgy与dat雷达数据样本及exe可执行程序，构成从源码到运行结果的完整工程。已有2015人学习下载，读者可据此梳理探地雷达数据处理的模块划分与算法实现思路，对照样本数据验证均衡、增益与滤波效果，并借鉴其可视化交互设计，为自主开发或二次改造同类处理工具提供参考。
1. 探地雷达数据处理程序：从原始回波到可解释剖面的落地路径
拿到一份探地雷达原始数据，最让人头疼的往往不是采集，而是采集之后那一堆二进制文件怎么变成能看的剖面图。探地雷达数据处理程序这类资源，核心价值就在于把去直流、增益补偿、滤波、背景去除、偏移这一串操作串成可复现的流程，而不是每次换一组数据就从头调参。它适合做道路检测、隧道衬砌、管线探测、考古勘探的一线从业者，也适合刚接触 GPR 数据、需要一套能跑通的最小工作流的新手。下面按「这份程序能干什么、怎么跑、参数怎么设、哪里容易翻车」的顺序拆开讲，中间会给出可直接抄的代码段和参数表，最后收在一个我常用的验证技巧上。
2. 数据读入与格式解析：把二进制头文件读明白
2.1 为什么格式解析是整条链路的地基
探地雷达数据常见的落盘格式有 GPR 厂商私有格式、SEG-Y、以及部分开源采集卡导出的自定义二进制。这类处理程序通常不会替你猜格式，而是要求你先明确三件事：道数、每道采样点数、采样时间窗。这三者决定了内存里那块二维数组的形状，一旦读错，后面所有滤波和增益都是在错误的数据上做无用功。我见过最常见的翻车，是把每道采样点数当成道数读，结果剖面出来是一条竖线，还以为是天线坏了。
程序里一般会有一个read_gpr之类的入口函数，负责按头文件里的字节偏移把数据块切出来。头文件字段的偏移量因厂商而异，所以资源里通常会附带一个字段映射表，你需要对照自己设备的头文件说明改偏移。下面这段是常见的读入骨架，假设头文件里已经解析出n_traces、n_samples、dt三个量。
import numpy as np def read_gpr(bin_path, n_traces, n_samples, dtype=np.int16, header_bytes=0): """按道-采样点二维结构读入原始 GPR 数据""" raw = np.fromfile(bin_path, dtype=dtype) # 跳过文件头，再按 (道数, 采样点数) 重塑 body = raw[header_bytes // np.dtype(dtype).itemsize:] data = body[: n_traces * n_samples].reshape(n_traces, n_samples) return data.astype(np.float32)逻辑说明：np.fromfile一次性把二进制读成扁平数组，header_bytes用来跳过文件头，reshape把一维流还原成「道 × 采样点」的矩阵。参数上，dtype必须和采集设备的量化位宽一致，16 位有符号是最常见的，若设备是 32 位浮点则要改成np.float32；header_bytes填错会让整幅剖面整体错位，表现为图像上下平移或出现周期性条纹。读入后建议立刻打印data.shape和data[:5, :5]，确认数值量级在合理范围，再往下走。
2.2 采样时间窗与深度换算的对应关系
读入之后紧接着要建立时间轴。探地雷达记录的是双程走时，深度换算依赖介质相对介电常数。程序里一般会提供一个time_axis和depth_axis的生成函数，前者用dt累加，后者用v = c / sqrt(eps_r)把时间换成深度。这里有个容易被忽略的点：dt的单位。很多头文件给的是纳秒，而程序内部按秒算，差三个数量级，深度轴会直接错到离谱。我一般会在读入后强制做一次单位归一，并在日志里打印dt的实际值，确认它在 0.05 到 2 纳秒这个常见区间内。
| 参数 | 常见取值 | 作用 | 设错的后果 |
|---|---|---|---|
| dt | 0.1–1.0 ns | 采样间隔 | 深度轴整体缩放错误 |
| eps_r | 4–9（土壤/混凝土） | 相对介电常数 | 深度换算偏差 20% 以上 |
| n_samples | 256–1024 | 每道采样点 | 剖面被截断或补零 |
| header_bytes | 依设备而定 | 文件头长度 | 数据整体错位 |
这张表建议直接贴在代码注释里，换设备时逐项核对。尤其是eps_r，很多人图省事填 6，结果在含水率高的路基上深度全错，后面解释层位时对不上钻孔，血泪经验就是这里来的。
3. 预处理链：去直流、增益与滤波的顺序问题
3.1 去直流与静校正为什么必须排在增益前面
原始道往往带一个直流偏置，表现为整道均值不为零，剖面上是一条亮带。去直流就是每道减去自己的均值，一行代码的事，但顺序很关键：如果先做增益再做去直流，增益会把直流偏置一起放大，去完直流后残余的偏置仍然可观。程序里通常把dewow（去低频）和dc_remove放在链路最前，我一般会先跑dc_remove，再跑一个短时窗的dewow，两者配合能把低频漂移压下去。
def dc_remove(data): """逐道去均值，消除直流偏置""" return data - data.mean(axis=1, keepdims=True) def dewow(data, window=15): """滑动均值去低频漂移""" kernel = np.ones(window) / window trend = np.apply_along_axis(lambda t: np.convolve(t, kernel, mode='same'), 1, data) return data - trend逻辑说明：dc_remove沿道方向求均值并减去，keepdims=True保证广播维度正确；dewow用滑动均值估计低频趋势再减掉，window控制去低频的尺度，太小会把有效信号一起削掉，太大则压不住漂移。参数上，window一般取采样点数的 1/20 到 1/10，256 点的话 15 到 25 比较稳。跑完这两步，剖面的整体亮度应该变得均匀，如果还是一条亮带，多半是window太小或者数据本身有饱和道，需要单独查那几道。
3.2 增益补偿：AGC 与手动增益怎么选
探地雷达信号随深度衰减很快，不做增益的话深部目标基本看不见。程序里一般同时提供 AGC（自动增益控制）和手动增益曲线两种。AGC 按滑动窗内的能量归一，适合快速出图看结构，但它会放大深部的噪声，做定量解释时不可靠；手动增益按深度分段给系数，可控但需要经验。我的习惯是先用 AGC 出一版预览，确认目标大致位置，再换手动增益做正式剖面。
def agc_gain(data, window=50): """滑动窗自动增益""" out = np.zeros_like(data) for i in range(data.shape[0]): trace = data[i] for j in range(len(trace)): lo, hi = max(0, j - window), min(len(trace), j + window) rms = np.sqrt(np.mean(trace[lo:hi] ** 2)) + 1e-6 out[i, j] = trace[j] / rms return out逻辑说明：对每个采样点，取前后window个点的均方根作为局部能量，用原始值除以它完成归一。window越大增益越平滑，越小则局部对比越强但容易出条纹。参数上，window一般取 30 到 80 个采样点，太小会在强反射处产生「光晕」，太大则深部增益不足。这段是纯 Python 循环，数据量大时会慢，实际用的时候可以换成向量化或分块处理，但逻辑不变。
3.3 带通滤波与背景去除的配合
滤波负责压掉带外噪声，背景去除负责消掉沿道方向一致的水平条纹（比如天线直达波和层间多次波）。两者顺序上，我一般先带通再背景去除，因为背景去除对高频噪声敏感，先滤波能让水平条纹估计得更干净。程序里通常用 Butterworth 带通，截止频率按天线主频来定，比如 400 MHz 天线，通带大概落在 200 到 600 MHz。
from scipy.signal import butter, filtfilt def bandpass(data, dt, low=200e6, high=600e6, order=4): """零相位带通滤波""" fs = 1.0 / dt b, a = butter(order, [low / (fs / 2), high / (fs / 2)], btype='band') return filtfilt(b, a, data, axis=1) def background_remove(data, method='mean'): """减去道间平均，压制水平条纹""" if method == 'mean': return data - data.mean(axis=0, keepdims=True) return data逻辑说明：butter生成带通系数，filtfilt做零相位滤波避免波形移位，axis=1表示沿时间轴滤。background_remove用所有道的平均作为背景估计再减掉，keepdims=True保证广播正确。参数上，low和high要落在天线带宽内，设得太窄会削掉有效反射，太宽则压不住噪声；order一般 4 阶够用，再高容易数值不稳定。背景去除对倾斜目标不友好，如果目标本身是斜层，减平均会把斜层也削掉一部分，这时候要改用滑动窗背景去除。
4. 偏移与成像：让绕射收敛到真实位置
4.1 为什么不做偏移剖面会「虚胖」
原始剖面里，点目标会表现为双曲线绕射，层位边界也会向外拖尾，看起来目标比实际大一圈。偏移的作用就是把绕射能量收敛回真实位置，让剖面接近真实几何。程序里一般提供 Kirchhoff 偏移或 f-k 偏移两种，前者实现直观、对速度模型要求相对宽松，后者快但对横向变速不友好。我一般先用常速 Kirchhoff 偏移出一版，速度取介质平均速度，看绕射是否收敛，再决定要不要上变速。
def kirchhoff_migrate(data, dt, dx, velocity, aperture=20): """常速 Kirchhoff 偏移，aperture 为偏移孔径（道数）""" n_traces, n_samples = data.shape out = np.zeros_like(data) for it in range(n_traces): for isamp in range(n_samples): amp = data[it, isamp] if amp == 0: continue t0 = isamp * dt for k in range(-aperture, aperture + 1): ix = it + k if ix < 0 or ix >= n_traces: continue # 双曲线走时：t = sqrt(t0^2 + (2*dx*k/velocity)^2) t = np.sqrt(t0 ** 2 + (2 * dx * k / velocity) ** 2) js = int(round(t / dt)) if 0 <= js < n_samples: out[ix, js] += amp return out逻辑说明：对每个输入样点，沿双曲线轨迹把能量累加到输出道对应时刻，aperture控制参与累加的道数范围。参数上，velocity是最关键的，取小了绕射收敛不足，取大了会过偏移出现「蝴蝶结」；aperture一般取目标深度的 1 到 2 倍对应的道数，太小能量收不全，太大引入噪声。这段是教学版实现，实际数据量大时要用索引预计算或 GPU 加速，但参数含义一致。
4.2 速度估计：别拍脑袋填一个数
偏移效果好不好，八成看速度。程序里一般会带一个速度扫描或双曲线拟合的小工具，用已知深度的目标（比如管线）反推速度。我的做法是找剖面上一个孤立强绕射，量它的顶点时刻和两侧同相轴的时差，用双曲线公式反算速度，再拿这个速度去偏移整幅剖面。如果整幅剖面速度变化大，就分段给速度，别一个数用到底。这里没有后悔药，速度错了偏移就是把能量搬错地方，越偏越乱。
5. 避坑与排查：五条踩过的坑
5.1 剖面整体发黑或发白
现象：出图后整幅剖面亮度异常，要么全黑要么全白。原因：增益系数过大或过小，或者去直流没做干净导致偏置被放大。解决：先关掉增益看原始数据量级，确认dc_remove生效，再把 AGC 的window调大一档，手动增益从 1.0 起步逐段加。
5.2 水平条纹去不掉
现象：背景去除后仍有明显水平亮带。原因：背景估计用了全道平均，而条纹只出现在部分道段，平均被稀释。解决：改用滑动窗背景去除，窗长取 20 到 50 道，或者先做道间中值滤波再减平均。
5.3 偏移后目标位置反而偏了
现象：偏移后绕射没收敛，或者收敛到错误位置。原因：速度填错，或者dx（道间距）单位与velocity不匹配。解决：核对dx是米还是厘米，速度用米/秒，两者统一后再跑；先用单目标反算速度验证，再全剖面偏移。
5.4 深部全是噪声
现象：增益后深部出现大片随机噪声。原因：AGC 在低能量区把噪声放大了。解决：换手动增益，深部系数封顶，或者先做带通滤波再增益，把带外噪声压掉再放大。
5.5 读入后道数对不上
现象：reshape报错或剖面道数明显偏少。原因：header_bytes偏移不对，或者dtype位宽与设备不符。解决：用十六进制工具看文件头前几十字节，对照设备说明确认偏移和位宽，先读前 10 道验证形状再全量读。
6. 验证与进阶：用已知目标反推参数是否可信
处理完一版剖面，怎么判断参数设对了？我常用的办法是找一个已知埋深的目标做闭环验证。比如某段道路已知有埋深 1.2 米的管线，就在剖面上量它的双程走时，用depth = v * t / 2反算速度，再拿这个速度回头检查偏移和深度轴是否自洽。如果反算速度和偏移用的速度差超过 15%，说明要么eps_r填错，要么时间轴单位有问题，整幅剖面的深度解释都不可信。
def validate_depth(two_way_ns, depth_m, c=0.3): """由已知埋深和走时反算相对介电常数""" v = 2 * depth_m / (two_way_ns * 1e-9) eps_r = (c / v) ** 2 return v, eps_r逻辑说明：two_way_ns是剖面上量到的双程走时（纳秒），depth_m是已知埋深，c取 0.3 米/纳秒。反算出的eps_r如果落在 4 到 9 之外，基本可以判定量测或单位有问题。参数上，c用 0.3 是为了和纳秒配套，换成米/秒时记得同步改two_way_ns的单位。
进阶一点的做法是把这套验证做成批处理：对每个已知目标跑一遍反算，统计eps_r的分布，取中位数作为全剖面的等效介电常数，再统一重算深度轴。这样比拍一个 6 靠谱得多。我现在的习惯是每换一组数据，先跑读入和去直流，打印形状和量级，再拿一个已知目标做速度闭环，确认无误才往下做增益和偏移。这套流程走顺了，换设备换工区也就是改几个头文件参数的事。希望帮到你。
本文还有配套的精品资源，点击获取
