简介：MatGpr2.0数据处理软件面向探地雷达数据分析的初学者与工程技术人员，基于MATLAB提供图形界面，用于完成GPR原始数据的导入、预处理、二维与三维成像、特征提取及结果解释，适合地质勘探、工程检测等场景下的教学与实操。压缩包共335个文件、约16.52MB，以174个m源码为主体，配有f90算法源码、dll与exe可执行组件、dzt/dat等示例数据，以及jpg/gif截图、html说明和pdf文档，便于对照界面与处理结果理解流程。已有850人学习，说明其在GPR数据处理入门、课程演示与案例复现方面具有一定参考价值。包内源码、示例数据与文档相互配合，读者可据此梳理噪声抑制、时深校正、二维与三维成像、目标解释及界面调用与参数设置的实现思路，并理解模块间的衔接关系，也可借助可执行程序快速跑通示例，对具备MATLAB基础者较易上手。

1. 探地雷达剖面的"第一公里"：MatGpr2.0 到底替你干了哪些活

刚拿到一段探地雷达原始剖面时，最费时间的不是调滤波参数，而是先判断这堆二进制里有没有可用的反射同相轴。MatGpr2.0 数据处理软件就是解决这个问题的：它把道头解析、零时校正、背景去除、增益、滤波、偏移、时深转换这些原本零散的操作，整合进一套 MATLAB 环境下的统一流程；包体里 MATGPRR_2 这一层是主流程入口，改参数、跑批、出图都在这一层完成。它适合手里已经有实测剖面、但被格式和参数卡住的人，也适合想把处理链固化成可复现脚本的人。下面按"先能读进来、再谈处理、最后谈出图"的顺序拆。

2. 环境与数据入口：MATGPRR_2 目录约定和首读校验的四个硬指标

探地雷达数据处理中有一个无法回避的现实：同一套流程，换一台机器，可能连第一步都跑不通。原因通常不是算法，而是路径里存在同名函数、MATLAB 版本不同导致图形句柄行为变化，或者数据读取时道头长度算错一位。MATGPRR_2 这一层之所以独立成一个目录，就是为了把主流程代码和数据、工具函数隔开，减少互相覆盖。下面先讲环境怎么搭，再讲数据怎么读进来，最后讲读进来之后如何用三条曲线判断这份数据是否值得继续处理。

2.1 MATLAB 版本与路径的最小可用配置

这类基于 MATLAB 的老工具包，最容易出问题的地方是界面回调。从 R2014b 之后图形句柄体系改过一次，凡是用了set(gcf,'...')老写法的地方，在新版本上可能报错或界面不刷新。我的习惯是先把主流程目录用genpath一次性加进来，然后立刻检查路径顺序，确认没有别的目录里有同名函数排在前面，导致调用到了错误的版本。

% env_check.m —— 最小可用环境自检 clear; clc; % 1) 记录当前 MATLAB 版本，出问题时第一时间能对上 v = ver('matlab'); fprintf('MATLAB 版本: %s\n', v.Release); % 2) 把主流程目录加到路径最前面 root = 'D:\gpr\MatGpr2.0'; addpath(genpath(fullfile(root,'MATGPRR_2'))); % 主流程入口 addpath(fullfile(root,'data')); % 原始数据目录 % 3) 检查路径顺序：MatGpr2.0 相关条目是否排在前面 p = strsplit(path, pathsep); idx = find(contains(p,'MatGpr2.0')); fprintf('MatGpr2.0 相关路径条目 %d 个，首个位于第 %d 位\n', numel(idx), idx(1)); % 4) 关键基础函数是否可见，避免被第三方工具箱顶掉 f = {'fft','hilbert','filtfilt','interp1','butter'}; for k = 1:numel(f) fprintf('%-10s -> %s\n', f{k}, which(f{k})); end

这段脚本的作用是把三个隐患提前明确处理。第一，版本信息要落下来，因为后面一旦出现"某张图显示不对"，第一件事就是确认是不是版本差异；第二，genpath会把子目录全加进去，但如果路径里已经存在旧版本的副本，MATLAB 是按顺序搜索的，先命中的赢，所以要看idx(1)是不是 1；第三，which的输出如果指向了某个不认识的工具箱目录，说明基础函数被覆盖了，这种情况在装了很多信号处理工具箱的机器上并不罕见，属于典型的隐蔽问题。

参数上只有root需要改，其余是自检输出。如果idx(1)不是 1，用rmpath把冲突目录移掉再跑一次。不要图省事直接addpath(genpath(root))加整个包，因为包里常有old、backup这类目录，里面放的是历史版本函数，一起加进去等于给自己留隐患。

2.2 把原始剖面读进来：道头、字长与道数对齐

探地雷达原始文件格式很多，常见的有 DZT、RD3、DT1、SEG2 这几类，它们的共同点是"道头 + 采样点"的定长结构，区别只在道头字节数和数据字长。读取最容易出错的地方不是文件打开，而是道数算错导致整幅剖面错位半道。我一般不用"文件大小除以单道字节数"直接取整就完事，而是把余数字节打出来看，余数不为零就说明道头长度或者字长猜错了。

% read_profile.m —— 定长道头读取 + 道数对齐校验 fn = 'line01.rd3'; ns = 512; % 每道采样点数 dt = 0.4; % 采样间隔，单位 ns hdr = 128; % 道头字节数 fmt = 'int16'; % 数据字长：int16 / int32 / single fid = fopen(fn,'rb','ieee-le'); % 小端序，多数仪器用这个 raw = fread(fid, inf, fmt); fclose(fid); per_trace = hdr/2 + ns; % 统一换算成 int16 个数 ntr = floor(numel(raw)/per_trace); raw = raw(1:ntr*per_trace); M = reshape(raw, per_trace, ntr).'; % 每行一道 data = M(:, hdr/2+1:end).'; % 去掉道头，转成 采样点×道数 fprintf('道数=%d 采样点=%d 冗余=%d\n', ntr, ns, numel(raw)-ntr*per_trace);

逻辑上有两步关键：hdr/2是把字节数换算成 int16 的个数，如果文件是 int32 或者 single，这里的换算系数要跟着变，不然整幅图会错位；reshape之后要转置，因为reshape是按列填充的，不转置会得到一道跨很多道的错误数据。最后的冗余输出是校验的关键：正常应该是 0，如果不是 0，说明hdr或fmt至少有一个不对，这时候不要强行floor之后就往下走，先把前 200 个 int16 打印出来看道头里有没有可识别的标记。

dt这个参数不能从文件里读到时，只能按仪器采集设置手填，它直接决定后面时深转换的深度标尺，填错 20% 的话，所有深度结论都要打折扣。这是我认为最该在第一步就确认死的参数。

2.3 首读校验：三条曲线判断数据能不能用

读进来之后，别急着去背景、加增益。先看三样东西：道间能量一致性、直达波位置漂移、噪声底幅度。这三样对应三个判断——数据有没有丢道、零时是不是稳定、后面做增益时该给多大动态范围。

道间能量一致性看的是每道能量沿道号的曲线。正常剖面这条曲线应该是平缓波动的，如果中间突然掉到接近零再弹回来，基本可以判定有丢道或者空道被读进来了。直达波位置漂移看的是每道最大绝对值出现的位置，探地雷达的地面反射一般在前几十个采样点，如果这个位置在道间来回跳几十个点，说明触发不稳定，零时校正必须逐道做而不是整体平移。噪声底看的是深部（比如后 20% 采样点）的标准差，这个数值决定了后面 SEC 增益指数的上限——噪声底高的时候把增益指数开大，只会把噪声一起放大。

% first_check.m —— 三条校验曲线 e = sum(data.^2, 1); % 每道能量 [~, p] = max(abs(data), [], 1); % 每道最大值位置 nb = std(data(round(0.8*ns):end, :), 0, 1); % 深部噪声底 figure('Position',[100 100 900 600]); subplot(3,1,1); plot(e); ylabel('道能量'); title('道间能量一致性'); subplot(3,1,2); plot(p); ylabel('采样点序号'); title('直达波位置漂移'); subplot(3,1,3); plot(nb); ylabel('标准差'); title('深部噪声底');

三条曲线出来之后怎么判断：能量曲线的极值比中位数低于 0.2 的道，标记为坏道，后面要么剔除要么插值补；直达波位置的标准差超过 3 个采样点，零时校正就得逐道算偏移量；噪声底如果沿道号有明显趋势（比如从前往后单调上升），说明耦合在变化，增益不能用一个全局指数。

3. 预处理链的先后顺序：直流偏置、零时校正、背景去除谁先谁后

预处理的三步——去直流、零时校正、背景去除——单独看都不难，难在顺序。顺序错了，后面每一步都在处理上一步留下的残留，结果就是参数怎么调都不对，最后怀疑是数据本身有问题。这一章把顺序为什么这么定讲清楚，再给一份可以直接照抄的实现。

3.1 顺序错了会怎样：三个可复现的反例

先说反例一：先做背景去除再去直流。背景去除用的是道间平均相减，如果每道都带一个恒定的直流偏置，这个偏置在道间平均里会被保留下来，减完之后直流不但没去掉，反而被"平均"成了一个更平滑的低频趋势，后面再想用逐道去均值清掉，就得多做一次。正确做法是先去直流，因为直流是逐道独立的、和邻道无关的，先把它清干净，后面的道间运算才不会被污染。

反例二：先做零时校正再做背景去除，但零时校正用的是整体平移而不是逐道对齐。整体平移之后，各道的直达波还是没对齐，道间平均相减的时候直达波没有被抵消，反而在剖面上留下一条随道号起伏的假同相轴。这条假轴看起来很"地质"，最容易骗人。

反例三：三者都做了，但背景去除的窗口开得太大。窗口拉到全剖面道数的一半，等于减掉了一个几乎不变的平均道，结果是整幅剖面的横向变化被抹平，弱反射和强反射一起被压。这类操作在实测数据上表现为"处理完看着干净了，但目标层不见了"，是比较典型的过度处理。

3.2 每一步的判据与写法

去直流的判据很简单：逐道均值应该接近零。做完之后重画一次mean(data,1)，如果还有明显的道间起伏，说明有低频漂移而不是纯直流，这时候要么加大高通滤波的截止频率，要么在去均值之前先做一次一阶差分。

零时校正分两种。触发稳定的数据，全剖面共用一个偏移量就够；触发不稳定的数据，必须逐道找起跳点。找起跳点不要用"最大绝对值位置"，那个位置在强弱反射变化时会跳，比较稳的做法是先做一次短时窗能量比，找到能量突增的第一个采样点。

背景去除本质是高通滤波在空间方向上的实现，滑动窗口的道数决定了被抑制的空间波长。窗口取 11 到 31 道是常见范围，具体看目标体的横向尺度——目标横向跨 5 道，窗口就不能小于 15 道，否则目标本身会被当成背景减掉。

% preprocess.m —— 去直流 -> 逐道零时校正 -> 滑动道平均背景去除 d = double(data); % (1) 逐道去均值：清掉直流偏置 d = d - mean(d,1); % (2) 逐道零时校正：以能量突增点为起跳点，向前补零对齐 win = 8; % 短时窗长度 en = filter(ones(win,1)/win, 1, d.^2); % 短时窗能量 [~, pick] = max(en(1:round(0.3*ns), :), [], 1); % 只在前 30% 找 ref = round(median(pick)); for k = 1:size(d,2) s = ref - pick(k); v = d(:,k); if s > 0 v = [zeros(s,1); v(1:end-s)]; % 向下搬，顶部补零 elseif s < 0 v = [v(1-s:end); zeros(-s,1)]; % 向上搬，底部补零 end d(:,k) = v; end % (3) 背景去除：沿道方向滑动平均后相减 w = 21; % 窗口道数，取奇数 bg = filter(ones(1,w)/w, 1, d, [], 2); d = d - bg;

这段代码有三处要留意。pick只在前 30% 采样点里找，是因为深部也可能有强反射，不限范围会挑错；补零而不是circshift，是因为循环搬移会把底部数据绕到顶部，制造出一条假同相轴；filter沿第 2 维做滑动平均，前后各有w-1个点的边界效应，如果目标体靠近剖面两端，那一段的背景去除是不完整的，要么裁掉，要么换用对称填充。

3.3 参数怎么定：窗口、阈值与迭代次数

窗口道数w是这三个参数里最需要试的。我的做法是先按目标体横向道数估一个下限，比如目标横跨 6 道，w至少取 18；然后从 21 开始往上加，每加 10 道看一次剖面，如果目标层的同相轴开始变淡，就退回上一档。这个试法比理论推导快得多，因为横向尺度在实测数据里本来就不好估。

短时窗长度win影响起跳点拾取的灵敏度。win太小，能量曲线毛刺多，pick会在相邻几个点之间来回跳；win太大，起跳点会被平滑到真实位置之后。8 到 16 个采样点是常见范围，对应dt为 0.4 ns 时大约是 3 到 6 ns 的窗长，和常见天线的脉冲宽度量级相当。

迭代次数指的是背景去除做几遍。做一遍是常规，做两遍在强背景干扰下有时候有效，但超过两遍基本就是过度处理了，弱反射会被一层层剥掉。判断依据是做完一遍之后道间平均的残余能量，如果残余已经低于噪声底，就没必要再来一遍。

4. 增益与滤波：弱反射信号从噪声里捞出来的参数边界

增益和滤波是整条处理链里参数最多、也最容易被滥用的环节。容易被滥用的原因是它们"立竿见影"——参数一调，图上立刻好看。但好看和可信是两件事，增益开过头会把噪声塑造成看起来像反射的结构，滤波带通开太窄会把真实信号的频谱削掉一半。这一章讲清楚什么时候该用哪种增益、滤波的边界在哪，最后给一份可以照着填的参数表。

4.1 SEC 与 AGC 的取舍：保振幅还是保可见度

SEC（球面与指数补偿）是按走时给一个解析增益，形式是t^alpha，它保留了道间振幅的相对关系，适合后面要做振幅分析、要做属性提取的场景。AGC（自动增益控制）是按滑动窗内的能量归一化，它把每段的振幅都拉到差不多的水平，视觉上最清楚，但道间振幅关系被破坏，做完 AGC 之后不要再用振幅去判断强弱反射。

选择依据很直接：如果目的是看结构、拾取同相轴、做解释，AGC 更省事；如果目的是做层位振幅对比、做衰减分析、算反射系数，只能用 SEC，而且alpha要尽量小。一个常见的折中是先用 SEC，把深部拉起来但不拉满，再在显示阶段用clim做百分位裁剪，把显示动态范围和数据处理分开。

alpha的取值和噪声底强相关。噪声底低的数据，alpha取 1.0 到 1.5 都能看；噪声底高的数据，alpha超过 1.5 之后深部基本全是放大后的噪声。判断方法是把增益后的深部区域单独画出来，如果出现大量横向不连续的高频条纹，说明增益过头了。

4.2 带通与 FK 的适用边界

带通滤波的前提是你知道天线的有效频带。常见做法是先把若干道的频谱平均一下，看主瓣在哪、-10 dB 带宽是多少，然后把通带设在主瓣两侧各留一点余量的位置。通带设得太窄，子波会变长，两个靠近的反射面会糊成一个；通带设得太宽，带外噪声进得来。4 阶 Butterworth 前后向滤波（filtfilt）是零相位的，不会引入时移，这是我默认会用的方式。

FK 滤波适合压制线性干扰，比如侧向反射、天线直达波的残留。它的边界在于：当目标体的同相轴斜率与干扰斜率接近时，FK 滤波会把目标一起压掉。所以做 FK 之前一定要先看一眼 f-k 谱，确认干扰的能量团和目标不在同一个区域，否则宁可不做。

% gain_filter.m —— SEC 增益 + 零相位带通 + 包络 t = (0:ns-1)'*dt; % 双程走时，单位 ns al = 1.2; % SEC 指数，先从 1.0 试起 g = (t + dt).^al; % 加 dt 防止 t=0 处除零 dg = d .* g; fs = 1/dt; % dt 单位 ns 时，fs 单位是 GHz [b,a] = butter(4, [0.15 0.55]/(fs/2), 'bandpass'); [sos,g0] = tf2sos(b,a); % 转二阶节，数值更稳 df = filtfilt(sos, g0, dg); % 零相位，不产生时移 env = abs(hilbert(df)); % 包络，用于看弱反射

这里有三个参数要说明。al是 SEC 指数，从 1.0 开始试是最稳的起点，噪声底低再加；[0.15 0.55]是归一化到 Nyquist 频率之后的通带边界，实际物理频率等于这个数乘以fs/2，比如dt=0.4 ns时fs=2.5 GHz，通带就是 0.375 到 1.375 GHz，要按自己天线的标称频带改；tf2sos那一步不是可有可无的，4 阶以上直接用[b,a]做filtfilt在低频端容易出现数值不稳定，转成二阶节级联之后就稳了。

包络那一步单独说一下。hilbert变换求包络在探地雷达里用得很多，它的好处是把振荡子波变成一个正的包络，弱反射在包络图上更容易看出来；代价是纵向分辨率下降，两个相距小于半个子波长度的反射面在包络上会合并。所以包络图适合做定性解释和层位追踪，不适合用来量厚度。

4.3 一条完整处理序列的参数表

把前面几章的操作串起来，一份可以直接照填的参数表长这样。这张表的价值不在于数值本身，而在于每一项后面写了"怎么改"——照抄只能跑通一次，知道怎么改才能用到自己的数据上。

表里最容易被忽略的是最后一行。介电常数取 9 还是取 16，深度结果差 33%，而很多人是直接用一个"常用值"填进去就出图了。有已知埋深的目标体（管道、电缆、已知界面）时，正确做法是用它反算eps_r，再拿这个值去标定整幅剖面，这比查表靠谱得多。

5. 常见问题排查：五类翻车现场的现象、原因与定位手法

处理链跑通不难，难的是结果不对的时候知道从哪儿查。下面五条都是我实际遇到过、并且能稳定复现的，按"现象 → 原因 → 解决"写，每条都给了定位手法，不是让你猜。

5.1 读取与坐标类问题

现象一：剖面整体倾斜，同相轴从左往右一路往下走，倾角还特别均匀。原因通常是道头长度猜错了，读取时把道头里的某个字段当成了采样点，导致每道都多读或少读固定的几个数，积累下来就成了均匀倾斜。解决方法是拿numel(raw)/per_trace的余数验证，余数不为零就说明长度不对；再打印前 128 个 int16，找有没有重复出现的固定标记（很多格式在道头里有同步字），用标记间距反推道头长度。定位手法：把第 1 道和第 50 道单独画出来叠在一起，如果两条曲线整体错开一个固定量，基本可以确认是这个原因。

现象二：整幅剖面沿道号方向有周期性的明暗条纹，周期正好等于某个整数道。原因是读取时把道头也当成数据算进了采样点里，每道末尾多出几个数，在图上表现为周期性的横向条纹。解决方法是回到hdr/2这个换算上，确认字节数和字长的对应关系没搞反。

5.2 处理链类问题

现象三：背景去除之后目标层不见了，但噪声也没了，图看着特别"干净"。原因是窗口w开得比目标横向尺度还小，目标被当成背景减掉了。解决方法是把w加到目标横向道数的 3 倍以上重新跑一遍，对比两次结果里目标同相轴的残余幅度。定位手法：把bg（减掉的那部分）单独画出来看，如果目标层在bg里清晰可见，那就实锤是窗口太小。

现象四：做完零时校正，直达波反而更乱了。原因是逐道拾取时用了最大绝对值位置，而深部强反射的幅度超过了直达波，拾取点跳到了深部。解决方法是限制拾取范围（只在前 30% 采样点找），并且改用短时窗能量比而不是单点绝对值。定位手法：把pick沿道号画出来，如果曲线中间有整段跳到 200 以后的，就是这个原因。

现象五：滤波之后剖面出现振铃，强反射两侧多出一串平行的弱同相轴。原因是带通通带设得太窄，滤波器阶数又高，时域响应变长。解决方法是把通带放宽，或者把阶数从 4 降到 2 再试。定位手法：拿单道做频谱，看主瓣宽度和通带的比值，比值小于 3 的时候就容易出振铃。

5.3 出图与导出类问题

现象六：同一份数据，imagesc显示正常，导出的图片却一片灰。原因是导出时用了print的默认渲染器，对某些插值方式支持不好。解决方法是显式指定渲染器，或者干脆用exportgraphics导出，并且导出前把clim固定死，避免导出时重新计算色标。定位手法：把caxis的返回值打出来，如果两次显示的值不一样，就是色标在自动重算。

还有一类容易被忽略的现象：剖面看着没问题，但一量深度就对不上。这多半不是处理的问题，而是dt或eps_r填错了，回到第 2 章和第 4 章的那两个参数上重新核一遍，比在滤波参数里反复试要快得多。

6. 从剖面到成果图：时深转换、包络出图与批量导出技巧

时深转换这一步，写法上有两套。一套是先算深度轴再显示，好处是坐标轴直接就是米，出图不用再改；另一套是显示时用走时轴，标注的时候再换算，好处是处理过程中不用反复改轴，但容易在最后出图时算错。我一般用第一套，因为成果图是给不熟悉走时的人看的，轴上直接是米最省沟通成本。

6.1 时深转换的两套写法与校准

% depth_convert.m —— 用介电常数把走时轴换成深度轴 eps_r = 9; % 干燥土体 4~9，湿土 9~16 c = 0.3; % 光速，单位 m/ns v = c/sqrt(eps_r); % 介质中速度，m/ns z = t*v/2; % 单程深度，m；t 是双程走时 figure('Position',[100 100 700 600]); imagesc([1 size(df,2)], z, df); colormap(gray); axis ij; xlabel('道号'); ylabel('深度 (m)'); clim(prctile(df(:), [2 98])); % 百分位裁剪，避免个别强点吃掉动态范围 colorbar;

clim那一步值得单独说。默认的色标范围取的是数据的最小最大值，探地雷达剖面里往往有一两个极强的直达波点，把色标拉到很宽，弱反射就全被压成灰色了。用 2% 和 98% 分位数裁剪，等于把最极端的 4% 数据截断，视觉动态范围立刻打开。这个操作只影响显示，不改数据，所以放心用。

eps_r的校准是这样做的：找一个已知埋深的目标，在剖面上读出它的双程走时t0，用eps_r = (c*t0/(2*h))^2反算，h是已知深度。算出来的值如果和查表差得远，先怀疑dt填错了，再怀疑目标深度记错了，最后才怀疑介电常数本身。

6.2 出图的三条显示建议

第一条，wiggle 图和灰度图叠着画。灰度图看整体结构，wiggle 看相位细节，两个一起画能同时照顾到"看趋势"和"看起跳"。叠的时候 wiggle 的横向缩放系数要按道间距调，缩放太大相邻道会互相压住。

第二条，包络图单独出一张。包络对弱反射的显示效果比原始剖面好很多，尤其是在深部衰减严重的情况下。但包络图的深度标注要和原始剖面一致，不要因为包络看起来"更平滑"就改了深度轴。

第三条，批量导出时把参数一起写进文件名。我吃过这个亏：一批剖面导出来，文件名只有线号，过两周再看完全想不起来当时用的是哪个alpha、哪个通带，只能重跑一遍。后来的做法是文件名里带上关键参数，比如line01_a1.2_bp015-055.png，一眼就能对上。

% batch_export.m —— 批量导出，参数写进文件名 al = 1.2; bp = [0.15 0.55]; tag = sprintf('a%.1f_bp%03d-%03d', al, round(bp(1)*100), round(bp(2)*100)); for i = 1:numel(lines) fn = sprintf('%s_%s.png', lines{i}, tag); exportgraphics(gcf, fullfile(outdir,fn), 'Resolution', 200); end

这段的关键是tag拼接，bp乘 100 再取整是为了避免小数点出现在文件名里，跨系统拷贝时少一类问题。Resolution给 200 是出图用的，如果要印刷再往上加，但注意别超过 600，再高对位图来说意义不大。

从一开始我把alpha和通带随手记在草稿纸上，到后来发现草稿纸永远找不到，再到现在的做法是参数必须进文件名，这个转变花了我不少时间。现在的习惯是：任何一批导出，先写tag，再跑循环，跑完对着文件名就能把当时的所有参数还原出来。希望帮到你。

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