探地雷达(Ground Penetrating Radar, 简称GPR)是一种利用高频电磁波(频率范围一般在1 MHz~4 GHz 之间)进行地下勘探的地球物理方法, 该方法通过发射天线向地下发射高频电磁波, 接收天线接收反射回地面的电磁波, 电磁波在地下介质中传播时遇到存在电性差异的界面时发生反射, 根据接收到电磁波的波形、振幅强度和时间的变化特征推断地下介质的空间位置、结构、形态和埋藏深度等信息^{[1, 2]}。低频探地雷达一般是指频率范围在1~50 MHz之间, 其探测深度较常规探地雷达有较大突破, 在一些特定环境下可以达到100~200 m以上^{[3, 4, 5, 6]}。而我国已经发现的陆域冻土区天然气水合物储层埋深多在200 m范围内^{[7, 8, 9, 10]}, 这使得应用低频探地雷达直接探测天然气水合物成为可能。为达到相对较深的探测深度, 低频探地雷达系统在设计时均采用较大的时间窗口(一般> 2 000 ns)及大尺寸天线及长收发距, 很难像高频GPR那样采用屏蔽装置压制噪声, 因此低频GPR信号在地下传播时受到能量更强的空气回波等噪声干扰, 深反射层回波较为微弱, 可以视为典型的非平稳、非线性信号。传统的傅里叶分析无法兼顾信号在时域和频域的不同特征, 所以诞生了Wigner-Ville分布、小波变换等改进的时频分析方法。但是小波函数需要针对所分析的信号选取合适的基函数, 才能得到较好的效果, 同时小波基受到测不准原理的限制, 无法在时域和频域同时为紧支撑区间, 易造成能量泄露。Wigner-Ville分布具有较好的视频聚焦特性, 但是对于多分量信号, 交叉项会产生虚假信号^{[11, 12, 13, 14]}。

对于非平稳、非线性信号, 较为直观的分析方法是使用具有局部特征的参量如瞬时频率、瞬时相位等。希尔伯特变换常常用来求解信号的瞬时频率和瞬时相位等特征, 但是很多信号经过希尔伯特变换后可能出现无法解释、缺乏实际物理意义的频率成分。Norden E. Huang等^{[15, 16, 17, 18, 19, 20]}在对瞬时频率进行深入研究后, 提出了一种新的时频分析方法— — 希尔伯特-黄变换(Hilbert-Huang Transform, HHT)。HHT包含两个关键步骤:经验模态分解(Empirical Mode Decomposition, EMD)和希尔伯特谱分析。EMD把信号分解为一系列本征模态函数(Intrinsic Mode Function, IMF)之和, 使得每个IMF的Hilbert变换得到的瞬时频率具有确定物理意义, 从而能够更加准确有效地分析信号的特征信息。在地球物理学领域, 基于EMD分解的Hilbert-Huang变换已经在地震信号分析、重力学、大地电磁测深法等领域得到越来越广泛的应用^{[21, 22, 23, 24]}。而在探地雷达信号处理领域, 已经有学者将HHT用于GPR剖面数据增强、GPR信号去噪、地下目标识别等方面^{[25, 26, 27, 28, 29]}, 这些研究都只是针对有屏蔽装置的高频GPR, 均没有涉及低频GPR及其信号处理等内容。

本文以HHT方法为基础, 研究低频GPR深部弱信号提取技术, 对低频GPR原始数据进行HHT变换, 可以得到各阶IMF其对应的瞬时属性信息, 从而有效区分地层反射信号和噪声, 在此基础上进行精细滤波, 压制深部噪声, 突出GPR深部弱反射信号。将该技术应用于陆域冻土区天然气水合物探测中, 能够有效分辨冻土层底界及其下分布的天然气水合物储层。

每个IMF刻画了信号在不同时频尺度上的特征, IMF需满足以下两个条件:①在整个信号中, 极值点和过零点的数量相等或最多相差一个; ②在任意点处, 由极大值点和极小值点定义的上下包络均值为零。任意信号s(t)可通过上述筛选分解为一系列IMF, 过程如下。

步骤一:将原始信号s(t)作为输入信号, 提取其所有局部极值点, 将全部极大值点和极小值点分别用样条曲线拟合出上包络和下包络, 使s(t)的所有点都包含在上下包络线之间。取上下包络的均值曲线 a(t), 从s(t)中减去a(t)可得到第一个被提取出的分量

e1(t)=s(t)-a(t)。(1)

步骤二:通常e₁(t)并不满足IMF条件, 需要把其作为输入信号重复步骤(a), 迭代k次后得到满足IMF条件的分量e_1k(t), 成为首个本征模态函数, 记为 f₁(t):

f1(t)=e1k=e1(k-1)-a1k(t)。(2)

步骤三:分解出第一个本征模态函数f₁(t)后, 从s(t)中减去f₁(t)得到余项r₁(t)

r1(t)=s(t)-f1(t)。(3)

步骤四:把r₁(t)作为新的输入信号重复步骤一~步骤三, 则可依次得到其余的IMF, 记为f₂(t), f₃(t), …, f_n(t), 余项记为r₂(t), r₃(t), …, r_n(t)。当余项r_n(t)变为单频函数或只含有一个极值点时, 上述分解过程停止。

这样, 经过EMD分解后, 原始信号s(t)被分解为若干IMF和余项r_n(t)的和:

这n个IMF可以视为原始信号s(t)的一组由信号本身定义的自适应基函数, 即不同信号经EMD分解会得到不同的基函数, 这与传统的傅里叶分解有着本质的不同。

GPR信号做上述EMD分解时, f₁(t)是最先被分解出的IMF, 包含了较多局部极值点, 因此反映了信号的高频特征, 噪声几乎都存在于f₁(t)中; 而随着尺度i的增加, f_i(t)最高频率越来越低, 而频率分辨率越来越高。从EMD分解过程可以看出, 信号和噪声在本征模态函数域中有不同的性态表现, 信号在IMF分量中通常随尺度增大而增大或不变, 而噪声在IMF分量中随尺度增大而迅速减小, 即噪声在各个尺度上的相关性较弱, 而信号在各个尺度上的相关性较强。因此, 可以使用IMF分量来构造IMF积检测器对分层信号进行检测, 通过多个尺度的IMF逐点乘积来凸显信号并抑制噪声。该方法仅仅需要信号脉宽这一先验知识, 可以对未知噪声分布、低信噪比情况下的GPR信号自适应检测。

对发射中心频率为15 MHz的GV6型号低频探地雷达实测GPR信号经过上述EMD分解得到的6个IMF如图1所示。从图中可以清楚地看到, IMF1是最先被分解出的IMF, 包含了较多局部极值点, 因此反映了信号的高频特征, 如噪声几乎都存在于IMF1中; 随着尺度i的增加, IMFi最高频率越来越低, 而频率分辨率越来越高。

得到了内在模式函数分量后, 对各阶IMF分量进行Hilbert变换:

y(t)=1πP∫-∞+∞c(τ)t-τdτ, (5)

其中, P为Cauchy主值, c(t)与y(t)组合成解析信号z(t):

图1

	Figure Option ViewDownloadNew Window
	图1 GPR信号EMD分解各阶IMF分量

z(t)=c(t)+iy(t); (6)

从而可以定义时变的幅值a(t)和相位θ (t):

a(t)=c(t)2+y(t)2, (7)θ(t)=arctany(t)c(t)。(8)

进一步可以定义瞬时频率:

ω(t)=dθ(t)dt。(9)

最终, 原始信号可以表达为

这里不考虑余量r_n, 因为它是一个单调函数或者一个常量。式(10)的实部定义为信号的Hilbert谱, 记为H(ω , t), 对其时间积分可得到信号的Hilbert边际谱:

h(ω)=∫0TH(ω, t)dt。(11)

HHT基于信号局部特征, 能自适应地分解得到平稳的IMF分量, 经Hilbert变换后得到的瞬时频率和Hilbert时频谱能够反映真实的物理过程, 可以很好地分析处理非线性、非平稳信号。图2为原始信号及图1各阶IMF分量对应的瞬时能谱图, 可以看到原始数据深部信号较弱, 经过HHT可以识别深部3 000 ns及3 800 ns左右两个高频震荡反射信号。

图2

	Figure Option ViewDownloadNew Window
	图2 各阶IMF分量对应的瞬时能谱

通过HHT方法对正演模拟数据处理与分析, 研究基于HHT方法的低频GPR弱信号提取技术。对于低频GPR而言, 能量较强的噪声干扰主要包括地表直达波、系统振铃以及地表反射体的反射回波^{[30, 31, 32, 33]}。直达波由于集中在探地雷达记录最初的很短时间段, 对我们识别深部地下介质反射影响不大, 地表反射体回波特征明显可以通过F-K滤波的方法加以滤除, 而系统振铃主要特征表现为近水平状贯穿整个剖面, 其频率特征复杂, 难以通过普通的频率域滤波技术滤除, 因此通过加入模拟系统振铃噪声的信号来研究HHT方法对低频GPR的效果, 模拟陆域冻土区天然气水合物成藏环境设计模型进行正演计算, 并在正演数据上加入系统振铃(System ring)噪声, 分析HHT方法的去噪效果。

设计模型为冻土层下砂岩层内含天然气水合物储层(图3a), 砂岩层上面覆盖有100 m厚度的冻土, 天然气水合物块体在砂岩内深度120~140 m, 水平方向110~138 m范围内, 冻土层介电常数为4, 电导率为0.01 mS/m; 砂岩层介电常数为6, 电导率为 0.1 mS/m; 天然气水合物层介电常数为4.5, 电导率为0.01 mS/m, 二维正演方法采用已经成熟的时域有限差分(FDTD)方法, 其正演结果如图3b所示。图4为加入系统振铃噪声的正演剖面, 可见受噪声影响GPR剖面的分辨率降低十分明显。

图3

	Figure Option ViewDownloadNew Window
	图3 正演模型及模拟结果

图4

	Figure Option ViewDownloadNew Window
	图4 加入噪声后的GPR正演模拟剖面

分别对加噪后的各道数据进行EMD分解, 得到各个IMF分量的剖面(图5)。低频GPR的EMD分解过程可以看做是一个频率分解过程, 整个剖面信号频率由高到低分解成各个IMF分量, 由于探测目标深度大, GPR信号在地下介质中传播时间长, 存在很严重的频散, 低频成分呈水平状干扰较多, 而GPR信号在冻土底界及水合物储层界面上的反射则呈现相对高频的特征, 因此IMF1剖面很好地反映了原始GPR剖面的反射信息。对比图3b的模拟结果, 发现其频率成分过高, 水合物层形态信息不够丰富; IMF2剖面包含了部分IMF剖面确实的信息, 也存在着大量低频干扰, 而IMF3~6剖面则涵盖了大部分引入噪声。因此, 将IMF1与IMF2剖面叠加以获得目标层位的大部分反射信息(图6a), 而其上的大部分干扰信息可以通过简单的频率滤波加以滤除, 最终得到还原度较高的图像信息(图6b)。利用HHT方法可以很好地区分系统振铃噪声与真实信号, 识别深部弱反射信号。3 冻土带天然气水合物探测应用

2016年7月在青海木里地区已经钻获天然气水合物样品的DK9井附近进行了低频GPR试验, 共布设两条剖面, 呈十字交叉状(图7), 其中Line1测线长600 m, DK9位于距离测线南端150 m处, Line2测线长600 m, DK9位于距离测线东端点170 m处, 采用的低频GPR仪器为中国科学院电子学研究所研制的伪随机超深GPR系统, 具体参数见表1。

这里重点介绍Line2线的处理方法及结果。首先进行单道分析。对DK9井处的单道数据进行HHT变换, 得到各阶IMF分量及其瞬时频谱(图8)。图8a中DATA曲线为原始信号波形, IMF1-5为各阶IMF分量曲线, R为余值。图8b中可以看到原始数据波形幅值强度变化不明显, 这是由于深部信号反射弱, 地质噪声较强所造成的。在IMF1和IMF2分量曲线上, 对天然气水合物储层均有不同程度的反映(IMF1更好)。从对应的瞬时频率分量分布上可以看出天然气水合物储层上方有明显的频率突变, 由高频振荡迅速降低至正常范围。通过对原始数据的频谱分析可以快速直观圈定大致的天然气水合物储层范围, 受噪声等影响需要进一步处理以便获得更为精确的储层信息。

图5

	Figure Option ViewDownloadNew Window
	图5 各阶IMF分量剖面及HHT结果

图6

	Figure Option ViewDownloadNew Window
	图6 HHT处理结果

图7

	Figure Option ViewDownloadNew Window
	图7 低频GPR试验剖面位置示意

表1

表1

表1 低频探地雷达主要参数 编码方式 8192组Golay编码 码元宽度 80 ns 每组相位编码长度 8个码元 发射峰值功率 1 kW 接收机灵敏度(含脉冲压缩) -130 dBm 混合采样间隔 5 ns 时变增益电路增益范围 -10~53 dB 线性平均次数 1/2/4/8/16 发射中心频率 12.5 MHz 脉冲重复频率 200 kHz 天线形式 12.5 m长电阻加载式偶极子天线 最大探测时窗 4 200 ns

表1 低频探地雷达主要参数

图8

	Figure Option ViewDownloadNew Window
	图8 DK9处单道GPR原始信号各阶IMF分量(a)及瞬时能谱图(b)

将IMF1与IMF2叠加后重构信号并对其进行中间滤波, 得到信噪比更高的单道信号, 并与DK9处测井电阻率曲线、岩心信息及水合物储层进行对比(图9), 均有较好的对应关系。根据录井资料, 确定约120 m处为永久冻土层底板, 从雷达图像可以看到该处信号存在明显的幅值和频率变化; 180~200 m之间层位为发现水合物实物的岩石层, 在GPR信号上也存在幅值和频率变化, 信号稍弱于冻土底板处且与其呈反相形态, 这一现象也符合理论计算的结果。以此处结果为标准, 追踪整个剖面, 得到二维GPR剖面解释结果(图10)。在二维GPR剖面上, 110~120 m范围内存在连续同相轴, 其下反射信号迅速减弱, 根据一维结果及正演模拟结果判断为冻土层底板, 在整个剖面上整体平缓, 局部稍有起伏。其下地层在探地雷达剖面上表现出复杂分层及起伏特征(图5红线以下灰色同相轴), 期间夹杂一些强反射信号区(图5黑色椭圆区域), 这些区域不受岩层岩性及分布影响, 呈现水平分布不连续特征, 这与近期地质研究发现的天然气水合物储层分别特征相吻合, 划为疑似含水合物层。

图9

	Figure Option ViewDownloadNew Window
	图9 DK9处单道GPR处理后信号与测井视电阻率及水合物层位对照

以上研究表明, 基于伪随机脉冲的低频大功率探地雷达在青藏高原冻土区探测探测深度可达200 m, 天然气水合物储层对于探地雷达信号有着明显幅值特征和频率突变特性, 可以作为冻土区天然气水合物储层天然气探地雷达识别标志。

图10

	Figure Option ViewDownloadNew Window
	图10 GPR剖面解释成果

1) 低频探地雷达采集时间窗口大, 深部信号受到地层衰减和频散的影响十分明显, 真实反射信号往往伴随大量噪声, 因此低频探地雷达信号可以视为典型的非稳态信号。引入对于非稳态信号处理效果非常好的Hilbert-Huang变换方法, 在其基础上研究了低频探地雷达深部微弱回波信号处理方法, 能有效地滤除信号中的噪声, 从而实现对深层微弱信号的检测和提取。

2) 对实测低频探地雷达数据的处理表明, 采用基于HHT处理方法的低频GPR探测技术可以用来直接探测冻土区天然气水合物储层以及冻土层底板。

(本文编辑:沈效群)