简介：一篇关于探地雷达图像数据处理及应用研究的PDF学术文献，面向地质探测、考古调查、道路质量检测等领域的科研人员与工程技术人员，旨在解决探地雷达信号受背景噪声干扰、目标识别精度不足等问题。资源为单个PDF文件，压缩包约335KB，内容精炼且结构完整，便于下载后直接阅读。文中系统分析了探地雷达单道数据的组成，将信号分为直达波、地表反射波、环境介质干扰、随机干扰和目标体反射波五类，并给出相应数据采集模型；在数据处理环节，重点阐述了利用均值法抑制背景干扰、运用HILBERT变换获取瞬时振幅、瞬时相位和瞬时频率特征图像的技术路径，同时涉及图像滤波、增强与分割等后续处理方法，并结合实际工程数据验证了有效性。目前已有346人学习，对于开展探地雷达信号处理、图像解释或相关算法研究的读者，是一份具有参考价值的技术文献。

1. 探地雷达图像数据处理：五种成份没分清，后面全是玄学

拿到一份探地雷达 B 扫描图像，最先看到的往往不是目标，而是一排排水平亮条纹——直达波、地表反射、天线耦合留下的振铃。这些固定干扰把双曲线状的目标反射信号压得只剩半截，靠肉眼很难断定底下到底有没有管、管线边界在哪。这篇论文 PDF 正是把这条链路讲清楚的东西：先拆单道数据成份，再用均值法把背景干扰减掉，最后用 HILBERT 变换取瞬时振幅、瞬时相位、瞬时频率三张剖面。适合做市政管线探测、隧道衬砌检测、道路结构层评估的人，也适合被“看图猜物”困扰的入门操作员。需要提醒的是，这份 PDF 是 2010 年的扫描版，OCR 经常把 HILBERT 打散成“HILBE RT”，不影响理解，但要有心理准备。

2. 背景干扰抑制的均值法：一条 axis 参数，决定目标是留下还是被抹掉

2.1 单道数据里到底有什么

论文第 1 章把单道 GPR 数据拆成了五种成份，这个拆法值得先抄下来：

列表如下：

• 直达波 a(t)：由发射天线直接进入接收天线，集中在记录最初的很短时间段，对识别深部介质影响不大。
• 地表反射波 b(t)：空气与地面阻抗突变产生，能量远大于地下回波，衰减慢，容易形成多次反射。
• 环境介质干扰 c(t)：高频成分，容易引起振铃效应。
• 随机干扰 r(t)：来自系统噪声和环境背景。
• 目标体反射波 s(t)：唯一想保留的部分。

加上采样离散化，最终数据就是 y(n) = a(n) + b(n) + c(n) + r(n) + s(n)。M 道采样点、N 道数据，组成一幅 M×N 的 B 扫描图像。

这个拆法的价值在于：a、b、c 三类在测线方向上走时几乎不变，属于固定背景；s 是双曲线形态，走时随道号变化。能不能干净地分离它们，决定了后面 HILBERT 变换出来的三张剖面是“特征增强”还是“噪声放大”。

2.2 均值法去背景：几行 Python 的事，但方向不能错

论文里对均值法的表述是：从每个 A 扫描中减去整个 B 扫描图像中所有相同双程走时的 A 扫描的平均。翻译成代码就是：

import numpy as np # data 形状约定为 (n_samples, n_traces) # 第0轴是双程走时采样点，第1轴是测线方向道号 n_samples, n_traces = data.shape # 沿测线方向（axis=1）对每个走时位置求平均 background = np.mean(data, axis=1, keepdims=True) # 逐道减去背景 cleaned = data - background

这段代码的关键是 axis=1。均值法假设固定干扰在测线方向上走时不变，所以“同一双程走时”的道数据取平均就能代表背景。而目标反射是双曲线，走时随道号变化，在固定走时上取平均时目标能量被摊到多个道里，均值很小，减完之后主体还在。

axis 选错是最常见的翻车点。如果数据读进来是 (n_traces, n_samples)，那 axis=1 减的就是时间方向，整个剖面都会被抹平。我一般拿到数据第一件事就是打印 data.shape，确认第 0 轴是时间采样。

2.3 均值法为什么有效，以及什么时候失效

均值法能成立，依赖一个前提：背景干扰在整条测线上近似平稳。直达波和地表反射波确实满足这个条件，天线耦合稳定的话，水平条纹在几百道里基本不变。这种情况下用整条测线平均，信噪比提升非常明显。

但测线一短就出问题。论文工程实例里测线约 4 m，目标埋深只有 20 cm，如果目标反射在较多道上都有明显响应，“平均背景”里就已经混进了目标的一部分。我在处理类似短测线数据时，通常改用滑动窗均值估计背景，只取当前道附近 5～11 道的平均：

from scipy.ndimage import convolve # 每 7 道一个滑动窗，沿测线方向平滑 kernel = np.ones((1, 7)) / 7.0 background_mv = convolve(data, kernel, mode='nearest') cleaned_mv = data - background_mv

为什么不用中值法？中值对少数异常道更鲁棒，但对弱双曲线目标同样不友好，目标如果占据的道数超过一半，中值背景也会把目标吃掉。均值法虽然简单，只要目标横向范围小于测线一半，结果一般都够用。

3. HILBERT 变换工程落地：瞬时振幅、瞬时相位、瞬时频率怎么算、怎么读

3.1 窄带信号假设与解析信号

论文把探地雷达信号当作窄带信号处理，这是 HILBERT 变换能落地的前提。窄带信号的瞬时频率有明确物理意义，宽带信号求出来的瞬时频率会在多个频率分量之间跳变，很难解释。

先回顾一下数学定义。实信号 f(t) 的 HILBERT 变换是 f(t) 与 1/(πt) 的卷积，频域响应是 -j·sgn(ω)。也就是说，变换后幅频特性不变，负频率成分做 +90° 相移，正频率成分做 -90° 相移。利用实信号与其 HILBERT 变换正交的特性，构造解析信号 z(t) = f(t) + j·f̂(t)，这个解析信号只包含正频率成分，且幅度是原信号正频分量的两倍。

实际用到的是三个量：

• 瞬时振幅 R(t) = sqrt(f²(t) + f̂²(t))，正比于该时刻雷达信号总能量的平方根。
• 瞬时相位 θ(t) = arctan(f̂(t)/f(t))，与反射波能量强弱无关。
• 瞬时频率 ω(t) = dθ(t)/dt，是瞬时相位的时间变化率。

这里插一句：GPR 脉冲本质不是严格窄带，但地下介质响应在局部是缓慢变化的，瞬时频率在目标边界以外还是稳定的。所以工程上仍然可以用，只是看到瞬时频率剖面里有零星跳变时，别急着怀疑代码。

3.2 scipy.signal.hilbert 的用法与参数

用 Python 复现这个流程非常简单，scipy 已经把 HILBERT 变换封装好了：

from scipy.signal import hilbert # 采样间隔：20 ns 时窗，512 个采样点 dt = 20e-9 / 512 # 沿时间轴（axis=0）对每道数据做 HILBERT 变换 analytic = hilbert(cleaned, axis=0) # 瞬时振幅：解析信号的模 inst_amp = np.abs(analytic) # 瞬时相位：解析信号的辐角，必须先做相位解缠 inst_phase = np.unwrap(np.angle(analytic), axis=0) # 瞬时频率：相位对时间的导数，注意 np.diff 会让道数少一行 inst_freq = np.diff(inst_phase, axis=0) / (2.0 * np.pi * dt) inst_freq = np.vstack([inst_freq, inst_freq[-1:, :]])

参数说明：dt 由论文的采样时窗 20 ns 和每道采样点数 512 计算得出，约 0.039 ns。axis=0 是因为瞬时振幅、瞬时相位、瞬时频率都是对时间定义的，必须沿时间轴做变换。np.unwrap 是相位解缠，把 [-π, π] 的跳变修正为连续相位，这一步漏掉的话瞬时频率剖面会布满飞刺。最后补一行是为了保持和原剖面行数一致，方便直接做图像对比。

3.3 三条瞬时剖面的物理含义与判读要点

瞬时振幅是反射强度的度量，空间分辨率更高，适合确定介质变化范围和目标体分布位置。论文图 4(c) 里，管线的分布范围就是在瞬时振幅剖面里看出来的。

瞬时相位反映时距剖面上同相轴的变化，因为与反射波能量强弱无关，所以弱反射也能显示出来，适合追踪地层变化和小断层。论文图 4(d) 里的“同相轴错乱”，就是新旧混凝土分界面在相位剖面上的表现。

瞬时频率反映介质岩性变化，对分界面更敏感。论文图 4(e) 里能看清双曲线顶部，便于确定埋深，同时看到混凝土与泥土的分界面。

三张剖面是从不同侧面看同一段数据，搭配着看才有意义。单独拿一张出来都容易误判：瞬时振幅只能告诉你“这里能量强”，瞬时相位告诉你“这里同相轴断了”，瞬时频率告诉你“这里介质变了”，三个信息对上，目标判定才站得住。

4. 200 MHz 管线探测复现：论文参数链与四条剖面的判读顺序

4.1 论文仪器参数速查表与设置逻辑

论文工程实例用的是意大利 IDS 公司的探地雷达系统，对龙阳路某处地下管线进行探测。参数链直接列出来：

200 MHz 天线的波长在空气中约 1.5 m，在地下按常见介电常数估算会缩短到 0.3～0.5 m，四分之一波长分辨率大约 10 cm 量级，和目标管径 12 cm 匹配。时窗 20 ns，按电磁波在地下约 0.1 m/ns 的传播速度估算，单程探测深度约 1 m，目标埋深 20 cm 在这个时窗的前三分之一段，反射信号完整，不会因为时窗太短被截断。

4.2 处理流程：先抑背景，再取瞬时信息

按论文的处理顺序，实际跑通的是四步：

第一步，读入原始 B 扫描数据，确认 shape 是 (n_samples, n_traces)。第二步，用均值法去掉背景干扰，得到 cleaned。第三步，对 cleaned 做 HILBERT 变换，得到瞬时振幅、瞬时相位、瞬时频率三张剖面。第四步，把原始剖面、去背景剖面、三张瞬时剖面并列排放，按顺序判读。

这一步一步走下来有一个好处：如果跳过第二步直接做 HILBERT，直达波和地表反射的能量远大于目标回波，瞬时振幅剖面里目标会被背景干扰淹没，瞬时相位也会被强反射的相位变化主导，三张瞬时剖面的价值就发挥不出来。背景抑制不是可选项，它是后面所有步骤的前置条件。

4.3 看图顺序：从双曲线到同相轴错乱再到分界面

论文图 4 的判读顺序值得记下来。原始剖面能看到管线的双曲线特征，但干扰存在时“无法准确断定”，这时候不要急着下结论。

去背景后的剖面里，目标反射信息得到加强，双曲线特征更清晰，但细节信号仍不完整。接着看瞬时振幅剖面，目标体分布范围一目了然；再看瞬时相位剖面，同相轴错乱标记出新旧混凝土分界面，反过来验证管线位置；最后看瞬时频率剖面，双曲线顶部更清晰，能确定埋深和混凝土与泥土的分界。

这个顺序的逻辑是从“有没有目标”到“目标在哪”再到“边界在哪”，一层层收窄。我实际处理数据时的习惯是：把五张图并排放在一个图像窗口里，鼠标从上往下扫，先扫原始剖面确认双曲线存在，再扫瞬时振幅圈定范围，最后用瞬时频率剖面量埋深。论文末尾列的 14 篇参考文献也值得顺藤摸瓜，其中复信号分析技术的几篇和这本 PDF 的处理思路一脉相承。

5. 避坑记录：处理 GPR 数据时四个容易翻车的环节

5.1 瞬时频率两端“发毛”：谱泄漏和端点效应

现象：处理完的瞬时频率剖面，图像最上面和最下面出现大量无规律的亮暗条纹，像头发丝一样乱跳，目标区域反而看不清。

原因：HILBERT 变换是全局变换，scipy 内部先做 FFT、把负频率置零、再 IFFT，对信号周期性很敏感。时间窗开头和结尾信号截断不连续，产生谱泄漏，解析信号在端点处失真，差分算子又把失真放大了。

解决：处理前对每道信号做 taper，把端点衰减到接近零：

from scipy.signal import windows taper = windows.hann(n_samples).reshape(-1, 1) cleaned_tapered = cleaned * taper

显示瞬时频率剖面时，再裁掉上下边缘各十几行，比硬着头皮看端点靠谱得多。

5.2 瞬时相位剖面出现水平横条纹：忘了相位解缠

现象：瞬时相位剖面里出现每道都有的水平带状跳变，看起来像地层分界面，打井验证却发现那里什么都没有。

原因：np.angle 或者 MATLAB 的 angle 函数返回的是 [-π, π] 的主值，相位实际值超过这个区间就会跳变回另一头，相邻采样点之间会差一个 2π。差分前不解缠，瞬时频率会出现巨大尖峰。

解决：对相位先做 unwrap 再做差分：

inst_phase = np.unwrap(np.angle(analytic), axis=0)

如果用的是 MATLAB，对应函数是 unwrap(phase, [], 2)，沿时间维解缠。这个错误最容易伪装成“发现新地层”的惊喜。

5.3 均值法把目标双曲线削没了：背景里混进了目标

现象：减完背景后，目标双曲线直接消失，只剩沿测线方向近似不变的残迹，看起来像是信号被整体减掉了。

原因：测线只有 4 m，目标反射在几十道里都有明显响应，整条测线平均的背景里已经包含了目标的一部分。目标越强、横向范围越大，被削得越严重。

解决：改用滑动窗均值，只用当前道附近几十道的窗口估计背景，或者先做增益均衡再减背景。窗口宽度一般取目标双曲线横向跨度的两倍以上，太窄会把双曲线当成背景减掉，太宽又回到了整条测线平均的问题。

5.4 频率剖面噪声大：带通滤波和空间平滑不能省

现象：瞬时频率剖面信噪比反而比原始剖面差，目标双曲线被噪声淹没，整张图都是雪花点。

原因：弱反射区域信噪比低，瞬时相位在无信号区是随机游走，求导后噪声被大幅放大。瞬时频率剖面本身就是对噪声最敏感的一张图。

解决：先做带通滤波再取瞬时参数，200 MHz 天线我一般取 50～600 MHz 通带，用零相位滤波避免相位失真：

from scipy.signal import butter, sosfiltfilt fs = 1.0 / dt # 采样率约 12.8 GHz sos = butter(4, [50e6, 600e6], btype='bandpass', fs=fs, output='sos') cleaned_f = sosfiltfilt(sos, cleaned, axis=0)

显示前对瞬时频率剖面做一次 3×3 中值滤波，也能把孤立噪声点抹掉，同时保住双曲线顶部的边缘。

6. 验证法：用合成双曲线模型把处理链路跑通，再上实测数据

6.1 合成双曲线模型：30 行代码跑通全流程

没有实测数据时，可以先合成一个点目标双曲线响应，把均值法和 HILBERT 变换全流程跑一遍，确认每张剖面输出符合预期，再拿这组参数去处理真实数据。这个方法我从那以后每次拿到新数据都强制走一遍。

import numpy as np from scipy.signal import hilbert # 与论文一致：20 ns 时窗，512 个采样点 n_samples, n_traces = 256, 64 dt = 20e-9 / n_samples # 目标参数 depth_idx = 80 # 双曲线顶点的采样点位置 x0 = 32 # 目标中心的道号 half_width = 30.0 # 双曲线展宽参数，越小弯曲越厉害 # 生成合成 B 扫描 data = np.zeros((n_samples, n_traces)) for xi in range(n_traces): peak = depth_idx + (xi - x0) ** 2 / half_width if 0 <= peak < n_samples: data[int(round(peak)), xi] = -1.0 # 负极性脉冲 data[20, :] += 0.6 # 模拟直达波固定背景 rng = np.random.default_rng(0) data += rng.normal(0, 0.05, size=data.shape) # 随机噪声 # 均值法去背景 background = np.mean(data, axis=1, keepdims=True) cleaned = data - background # HILBERT 变换 analytic = hilbert(cleaned, axis=0) inst_amp = np.abs(analytic) inst_phase = np.unwrap(np.angle(analytic), axis=0) inst_freq = np.diff(inst_phase, axis=0) / (2.0 * np.pi * dt) inst_freq = np.vstack([inst_freq, inst_freq[-1:, :]])

跑完后检查：瞬时振幅剖面在 depth_idx 附近应该有明显的局部峰值，且峰值横向展宽和输入的 half_width 一致；瞬时相位剖面能看到双曲线同相轴连续变化；瞬时频率剖面在顶点附近稳定，边缘有少量跳变但不影响整体判读。如果这三条都满足，说明处理链路是通的。

6.2 两条快速自检的习惯

第一个习惯是算信噪比。选目标区和背景区，分别求能量比：

# 目标区：顶点附近 40 个采样点 × 16 道 target_region = cleaned[depth_idx-20:depth_idx+20, x0-8:x0+8] # 背景区：浅层直达波之后的安静区 bkg_region = cleaned[30:60, :] snr = 10 * np.log10(np.mean(target_region**2) / (np.mean(bkg_region**2) + 1e-12)) print(f"去背景后目标区相对背景信噪比: {snr:.2f} dB")

第二个习惯是永远保留原始剖面。所有处理都从原始数据派生，不在原数组上原地修改。这样一旦瞬时剖面出现无法解释的异常，随时可以回到原始剖面排查是处理逻辑的问题还是真实地质异常。

真正让我把这个流程固定成习惯的一次经历，是在瞬时相位剖面里把没解缠的相位跳变认成了混凝土分层，先入为主以为管线下面还有一层空洞，后来在合成模型上一跑才发现是 unwrap 漏了。从那以后，我每次处理 GPR 数据都先跑一遍合成模型，确认端点、相位解缠、背景抑制都正常，再切到实测数据。这个方法不复杂，但能帮你少走大半截弯路，希望帮到你。

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