ISSN 
下载:
-
井中雷达系统尚未有统一的国际化标准[1-4]，井中雷达数据传输速度有限且数据很容易因为复杂的环境而被“污染”，如何从原始数据中提取出有效数据，并且对井中仪器进行实时控制及监控[5-7]，是井中雷达系统数据处理所面临的首要问题。
在信号预处理方面，井中雷达信号因为工作环境复杂，遇到的干扰十分剧烈，同时由于其非线性、非稳态、超宽带的特点，导致其使用傅里叶变换、小波变换等常规分析方法难以实现信号分析[8-10]。井中雷达信号预处理工作需要将雷达对地下异常目标体的回波信号进行提取分析，由于目标反射较弱，或是探测深度较深导致回波信号衰减较大，导致信号强度很弱，同时又由于其反射环境复杂，容易受系统、井壁、井液等影响，导致杂波干扰非常剧烈[11-12]。因此，如何将复杂的回波数据分离出噪声和有用信号，获取真实雷达反射回波，是信号预处理甚至雷达硬件系统设计工作所面临的最大问题，具有重要的研究意义。
文献[13]将低维离散余弦变换(discrete cosine transform, DCT)特征向量与支持向量机(support vector machines, SVM)分类方案结合使用，采用Wiener滤波技术对探地雷达信号进行了自适应预处理，研究结果表明该方法能有效增强信号的分类精度，在计算效率和准确性方面能取得良好的性能。文献[14]发现探地雷达信号中有较强的散射波以及由色散引起的传播噪声，因此采用二维频域响应滤波器和基于传播常数的滤波来抑制其散射，就信号幅度水平而言，该滤波器在特定应用下信噪比有一定提高。文献[15]提出了一种反滤波方法用于消除探地雷达信号的子波衰减和频散，该方法认为地层介质的边界条件系数是随机数，充分利用了最小相位滤波器作为其等效滤波器，有效提高了探地雷达的探测深度，并且能够进一步提高探测分辨率。
另一方面，通过EMD (empirical mode decomposition method)对信号进行模态分解，这种分解方式不同于其他滤波和小波变换等方法，EMD分解对信号具有自适应性，它的优点是不会运用任何已经定义好的函数作为基底，而是根据所分析的信号自适应的生成固有模态函数[16]。文献[17-18]先后提出了EEMD (ensemble empirical mode decomposition)算法、高维集合EMD分解(modified ensemble EMD, MEEMD)，通过向原始数据中加入噪声，能够有效减小经验模态分解中存在的混叠问题，同时也将其推广到高维处理。文献[19-25]充分利用了模态经验分解的特性，根据EMD分解进行了多分辨率以及时频分析方法的研究，进一步拓展了EMD分解技术在探地雷达信号处理中的降噪方法和应用前景。
以上方法仅运用到探地雷达的信号处理中，并没有应用到井中雷达。本文基于EMD分解的信号预处理方法，对希尔伯特−黄变换算法进行了阐述，通过构造函数的方式验证算法，并对算法进行了井中雷达的信号预处理。
-
井中雷达的基本工作原理如图1所示。其接收信号可看作由发射信号与系统各个模块、反射信号以及周围环境等作用产生的，通过式(1)进行线性模型描述：
$$ w(t) = d(t) + c(t) + s(t) + n(t) $$ (1) $$ d(t) = {f_s}(t) * {f_{a1}}(t) * {f_c}(t) * {f_{a2}}(t) $$ (2) 式中，
$ d(t) $ 表示发射和接收天线间的直达波分量；$ {f_s}(t) $ 是经过脉冲源产生的信号在天线上产生的激励信号，会由于信号拖尾、波形畸变等对回波信号产生一定的影响；$ {f_{a1}}(t) $ 是该脉冲信号经过发射天线时所产生的冲激响应；$ {f_{a2}}(t) $ 是接收端接收天线对直达波信号的冲激响应；$ {f_c}(t) $ 是由于收发天线间的互耦效应所产生的互耦响应，主要包含了天线互耦作用$ {f_{c1}}(t) $ 以及井内泥浆和雷达仪器周围冲洗带影响$ {f_{c2}}(t) $ 。一般而言，井中雷达在实际工作中，雷达收发天线间的相对位置固定，直达波仅与仪器周围、收发天线之间的介质情况有关。因此，在探测过程中，直达波分量仅受井中填充井液影响，在一定井深范围内直达波分量是相对稳定的。
图 1 井中雷达模型
c(t)主要是电路耦合分量，仅针对收发一体的井中雷达仪器，由于收发一体的井中雷达接收天线与发射天线合而为一，通过 TR 组件对天线收发状态的切换，导致源与采集之间隔离较为困难使电路耦合分量，随电路工作稳定性呈一定波动。
$ s(t) $ 为雷达波信号经过地层传输，遇到目标反射而被接收端获取的信号回波分量，也是井中雷达探测所感兴趣的信号回波分量：$$ s(t) = {f_s}(t) * {f_{a1}}(t) * {f_g}(t) * {f_t}(t) * {f_{a2}}(t) $$ (3) 式中，
$ {f_g}(t) $ 是反射、折射、散射所引起的地层冲激响应；$ {f_t}(t) $ 为探测地层中目标的冲击响应，一般是由地层中孔洞、裂缝、埋藏物等产生，受目标的大小、形状、电磁特性以及传播地层的电磁特性影响，从而表现出不稳定的波形，产生目标物的相应波形，是井中雷达探测的目标分量，而雷达信号预处理所要达到的目的便是获取该信号的所有分量。n(t)为加性噪声分量，其来源较为复杂，而因为井中雷达的特殊性，其主要来源是系统热噪声。不同介质由于电磁特性不同，所得到的回波信号也不同。井中雷达回波信号受到的干扰来源复杂，极易受到井中及井周复杂的介质环境影响，且由于井中雷达器件的局限性，其接收到的信号信噪比较低，受到井周地层介质的干扰强烈，波形畸变更为严重。因此，通过合适且针对性较强的信号预处理手段来进行波形修正，降低环境、器件干扰，提高系统信噪比，进行井中雷达信号预处理从而获得更好的雷达探测结果。
-
井中雷达数据由3种形式构成：A-扫、B-扫、C-扫，分别对应雷达探测一维、二维、三维数据，表现形式分别是回波波形、波形矩阵灰度图、灰度图阵列。
井中雷达A扫数据由一次回波采样构成，能够描述在某一井深处，垂直井向的特定水平方向内的回波信号。其数学表达式使用柱坐标系描述为：
$$ f(r) = u( {{r_i},{\varphi _j},{z_k}} )\begin{array}{*{20}{c}} {}&{} \end{array}i = 1,2, \cdots ,N $$ (4) 式中，j、k为固定值，指示该回波信号的方向角与深度；i为时间采样点数；N为时窗长度。
B扫数据是由一定井段范围内的多个A扫数据组成，描述的是在井段范围内，某一特定水平方向内的回波信号。其数学表达式使用柱坐标系描述为：
$$ f(r) = u( {{r_i},{\varphi _j},{z_k}} )\begin{array}{*{20}{c}} {}&{i = 1,2, \cdots ,N\;\;\;\,k = 1,2, \cdots ,M} \end{array} $$ (5) 式中，j为固定值，指示特定水平方向；k为深度采样点数；M为某深度范围内采样波数。
C扫数据是井中雷达多个通道、不同方向，天线在相同井深处，多个方向内的回波信号，由雷达仪器在固定深度处由固定方向的天线旋转采样多次回波信号获取，或者由多个固定方向接收天线同时接收获得。其数学表达式使用柱坐标系描述为：
$$ \begin{split} &f(r) = u( {{r_i},{\varphi _j},{z_k}} )\;\;\;\;i = 1,2, \cdots ,N\\ &k = 1,2, \cdots ,M\;\;\;\;\;\;\;j = 1,2, \cdots ,P \end{split} $$ (6) 式中，j为雷达通道编号；P为雷达仪器通道数。
-
井中雷达的数据后期处理流程主要分为数据预处理、信号预处理、成像与反演3个阶段。其中，信号预处理主要是将信号进行去噪、滤波、去直流分量、增益控制处理，目标是剔除所有干扰信号，仅留下有用的回波信号。常规的方法有平均滤波、带通滤波、滑动窗口滤波、自动增益控制、小波变换等。
井中雷达预信号处理过程中，滑动平均滤波有时能更好地去噪，有效抑制周期性和渐变状的干扰，如由系统温度引起的采样时钟偏移，但滑动平均滤波不适用于由脉冲干扰所造成的噪声，即由时域波形的脉冲采样值误差，如采样噪点、环境干扰等，因此不适用于随机强噪点干扰严重的场合。同时，与去背景噪声算法一样，对平行于井向的平面或线性信号是无能为力的，需要针对不同的场景来应用不同的处理方法。一般而言，井中雷达信号在传播过程中，几乎不会产生频率偏移，因此，根据雷达的发射源的频率情况，对接收数据进行相应的带通滤波，滤波系数主要取决于发射信号的频谱。
但自然信号通常都是非线性和非平稳的，因此，在处理这些非线性和非平稳时域信号时，无论是傅里叶变换还是小波变换都会对信号引入线性分解，不是最佳方法。相较而言，对于井中雷达信号而言，由于信号的强非线性以及非平稳性，EMD无疑是非常有效的分析方法。通过对信号的分解及时频域分析，对信号根据其信号自身的特点，自适应地分解信号，为信号的时频域分析提供核心思路。
-
希尔伯特−黄变换算法是目前热门的信号分析算法，属于一种经验方法，其主要过程分为模态经验分解EMD算法和希尔伯特谱分析两个部分，通过对信号进行时域分解，从时频域对信号进行详尽的解析。
-
EMD的基本思想主要在于对所有复杂信号而言，其时间序列都可以分解成具有不同特征尺度的有限个本征模函数(intrinsic mode function, IMF)，即IMF分量。每一个IMF具有局部正交、自适应的特点，用于分析非平稳信号时，有效避免了不同分量之间随时间变化的相同频率部分互相干扰，大大增加了非平稳信号的分析预测和解释能力。
EMD的原理来源于一个简单的假设，即任何信号都由不同的IMF组成，每个IMF在一个分离的时间尺度上代表一个嵌入的特征振荡。EMD分解具有的局部自适应性，适用于非线性以及非平稳信号的分析。EMD的出发点是考虑局部振荡水平上的振荡信号，并对作为新信号的慢振荡分量进行迭代，如式(7)所示：
$$ {{s}}({\text{t}}) = {s_{{\rm{fast}}}}(t) + {s_{{\rm{slow}}(t)}} $$ (7) 计算流程如图2所示。
该方法旨在对非平稳和非线性信号进行分解，以展示其隐含的准周期性和特征。它是一种局部自适应方法，以数据作为驱动，因此更适合于非线性和非平稳数据的分析。
图 2 EMD分解迭代流程图
常用的数据分析方法中，大部分基于傅里叶变换进行，但用傅里叶变换进行数据分析时，通常需要满足两个条件：信号必须是平稳的，并且在给定的时间跨度内是线性的，然而一般来说，物理信号不满足这些条件。EMD分解方法将信号分解为振幅和频率调制函数，也就是固有模态函数IMF。EMD分解可以在自适应时间−频率−幅度空间中处理非线性和非平稳信号的数据，将信号进行分解为IMF和残差，给定方程如下：
$$ I(n)=\sum _{m=1}^{M}{\text{IMF}}_{m}(n)+{\text{Res}}_{M}(n) $$ (8) 式中，
$ I(n) $ 是多分量的待分解信号；$ {\text{IM}}{{\text{F}}_m}(n) $ 是第m个固有模态函数；$ {\text{Re}}{{\text{s}}_M}(n) $ 为对应于M个IMF的残差，或者称为趋势(mean trend)。EMD算法计算过程概述如下，以s(t)作为输入信号：
1) 找出s(t)局部极大值和极小值点，获得极大值点和极小值点的序列，分别将极大值点和极小值点进行3次样条插值，此时，获得s(t)的上包络线和下包络线，两个包络线互相不对称；
2) 对极大值包络线和极小值包络线取平均值，此时获得了包络均值，将原始信号减去均值包络线，此时获得第一个分量：
$$ {h_1}(t) = s(t) - {m_1}(t) $$ (9) 3)将
$ {h_1}(t) $ 作为输入重复进行步骤1)、步骤2)获得$ {h_2}(t) $ ，依此类推；4)重复上述步骤，使用相邻两次分解分量的标准差作为停止准则，一般取值为
$ 0.2 \leqslant S \leqslant 0.3 $ ，S计算如下：$$ S=\sum _{t=0}^{T}\frac{{\left|{h}_{k-1}(t)-{h}_{k}(t)\right|}^{2}}{({h}_{k-1}{(t)})^{2}} $$ (10) 此时获得
$ {h_p}(t) $ ，当再次重复步骤1)、2)时，将会获得均值包络为$ x = 0 $ 的一条直线或接近零的直线。此时，$ {h_p}(t) $ 的上下包络线具有良好的对称性，继续筛选下去，信号不会有太多的变化，因此停止筛选。5)当满足停止准则时，
$ {h_p}(t) $ 满足IMF分解条件，此时第一次EMD分解过程完成，$ {h_p}(t) $ 为第一个IMF，记为$ {c_1}(t) $ 。将原始输入信号减去$ {c_1}(t) $ ，获得剩余量$ {r_1}(t) $ ：$$ {r_1}(t) = s(t) - {c_1}(t) $$ (11) 6)将
$ {r_1}(t) $ 作为新的信号输入，重新执行步骤1)～步骤5)，获得新的余量$ {r_2}(t) $ 以及第二个模态$ {c_1}(t) $ ，如此重复n次，依次获得$ {c_3}(t),{c_4}(t), \cdots ,{c_n}(t) $ ，此时剩余的余量$ {r_n}(t) $ 成为单调函数，无法再进行IMF分解，整个EMD分解过程完成。由上述步骤可得：
$$ s(t)=\sum _{k=1}^{n}{c}_{k}(t)+{r}_{n}({t}) $$ (12) 此时，原始信号被分解成为n个模态分量
$ {c_k}(t) $ 以及剩余一个单调趋势函数分量$ {r_n}(t) $ 。至此，模态经验分解完成。 -
通过EMD分解后，将原信号分解为一组单分量信号IMF的组合，同时对算法进行模拟信号的分解验证，接下来对其分解结果进行进一步的希尔伯特谱分析。
将各IMF分量分别进行希尔伯特变换，记为
$ {\boldsymbol{H}}_{{c_k}} $ ，此时就得到了各个IMF的瞬时特征量，因此原信号可以表示为：$$ s(t)=\sum _{k=1}^{n}{c}_{k}(t){\text{e}}^{\text{j}\int {w}_{i}(t)\text{d}t} $$ (13) $$ s(t)=\sum _{k=1}^{n}{c}_{k}(t){{\rm{e}}}^{{\rm{j}}\int {w}_{i}(t)\text{d}t} $$ (14) 也即希尔伯特谱，可以看出，希尔伯特谱与傅里叶谱相比，傅里叶谱的推广更具有一般化的特点，希尔伯特谱更具有频率的瞬时特性。
此时对每一个IMF分量进行解析信号的构造，构造函数为：
$$ {Z_i}(t) = {c_k}(t) + {\text{j}}{H_{{c_k}}} $$ (15) 再使用欧拉公式对构造函数进行解析：
$$ {Z_i}(t) = {A_i}(t){{\text{e}}^{\text{j}}}\theta (t) $$ (16) 式中，
$ {A_i}(t) $ 即为分量信号的幅度序列，而通过计算能得到每个IMF的瞬时频率：$$ {W_i}(t) = \frac{{{\text{d}}\theta (t)}}{{{\text{d}}t}} $$ (17) 最终，通过EMD分解的解析函数的构造提取，即可绘制出原始信号的“时间−频率”分布，并且在每个时间点上，都能找到每个平稳态分量的频率值。
-
EMD分解的二维表现形式为二维模态经验分解(BEMD)，也叫2D EMD。BEMD同样需要寻找B扫数据的局部极大值点以及极小值点，分别对极大值点和极小值点进行后续插值，获得极大值曲面和极小值曲面，而提取每一个BIMF也需要多次迭代，其计算流程与一维EMD分解算法大致相同。
由于井中雷达特性所造成的非线性、不平稳、频带广的干扰噪声，对于B扫数据使用局部适应的二维EMD分解方法，可以有效去除非线性干扰，进一步保证数据质量。二维EMD分解主要用于雷达B扫数据滤波去噪，通过对雷达B扫数据进行二维EMD分解得到不同尺度特征的BIMF分量。由于在提取过程中引入了极值点距离，从而定义了局部尺度，同时在分量提取时，由于筛选条件的限定，由BEMD方法获取的图像，分解成为一组由“精细”至“粗糙”的图像组。首先，分解出的是比较细节的部分，包含了较为高频的分量图像信息，之后分解出的图像逐渐尺度变大，显示了较为低频的分量图像信息，最终剩余部分是最大尺度的趋势分量图像信息。BEMD方法对雷达B扫数据具有一定的自适应性，计算过程大致与一维相似。
二维EMD分解涉及到的插值方式与一维较为不同，本文使用的二维插值方式为多重二次曲面(multi-quadric, MQ)插值，MQ插值方案非常适用于二维包络插值，该插值方式主要用于对地理表面以及重力和磁场表面进行插值，已被广泛应用于地球科学和其他领域。MQ方法是一种使用径向基函数(radial basis function, RBF)的散点数据全局插值方案。对于散点数据
$ ({x_i},{y_i}) $ 极值，MQ方法使用方程式：$$ H(X)=p(X)+\sum _{j=1}^{n}{\lambda }_{j}\phi ({X}_{j})\begin{array}{cc}& j=1,2,\cdots ,n\end{array} $$ (18) 式中，H(X)是极值点
$ ({x_i},{y_i}) $ 的值；p(X)为低价多项式，$ {{p}}({{X}}) = {c_0} + {c_1}x + {c_2}y $ ；$ {\lambda _j} $ 是待定参数；$\phi = [{({X_i} - {X_j})^2} + {R^2}]^{1/2}$ 为内插常数。将式(18)在矩阵中表示，可以写成：
$$ \left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{\boldsymbol{P}} \\ {{{\boldsymbol{P}}^{\text{T}}}}&{\bf{0}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {\boldsymbol{\lambda}} \\ {\boldsymbol{C}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\boldsymbol{H}} \\ {\bf{0}} \end{array}} \right] $$ (19) 其中：
$$ \begin{split} & \quad {{\boldsymbol{A}} = \phi \left| {{X_i} - {X_j}} \right| = } \\ & {\left[ {\begin{array}{*{20}{c}} {{d_{11}}}&{{d_{12}}}& \cdots &{{d_{1n}}} \\ {{d_{21}}}&{{d_{22}}}& \cdots &{{d_{2n}}} \\ \vdots &\vdots & \ddots & \vdots \\ {{d_{n1}}}&{{d_{n2}}}& \cdots &{{d_{nn}}} \end{array}} \right]} \end{split} $$ (20) 式中，
$\phi $ 为内插常数，矩阵中的每一个元素为：$$ {d_{ij}} = {\left[ {{{\left( {{x_i} - {x_j}} \right)}^2} + {{\left( {{y_i} - {y_j}} \right)}^2} + {R^2}} \right]^{\frac{1}{2}}} $$ (21) 元素P为：
$$ {\boldsymbol{P}} = \left[ {\begin{array}{*{20}{c}} 1&{{x_1}}&{{y_1}} \\ 1&{{x_2}}&{{y_2}} \\ \vdots & \vdots & \vdots \\ 1&{{x_n}}&{{y_n}} \end{array}} \right] $$ (22) $ {\boldsymbol{\lambda}} $ 、$ {\boldsymbol{C}} $ 、$ {\boldsymbol{H}} $ 分别为：$$ {\boldsymbol{\lambda}} = {[{\lambda _1},{\lambda _2}, \cdots ,{\lambda _n}]^{\text{T}}} $$ (23) $$ {\boldsymbol{C}} = {[{c_0},{c_1},{c_2}]^{\text{T}}} $$ (24) $$ {\boldsymbol{H}} = {[{H_1},{H_2}, \cdots ,{H_n}]^{\text{T}}} $$ (25) 在式(19)中，只有矩阵
$ {\boldsymbol{\lambda}} $ 和C是未知的，并且$\left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{\boldsymbol{P}} \\ {{{\boldsymbol{P}}^{\text{T}}}}&{\bf{0}} \end{array}} \right]$ 是实对称矩阵，因此，参数的解可以表示为：$$ \left[ {\begin{array}{*{20}{l}} {\boldsymbol{\lambda}} \\ {\boldsymbol{C}} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{\boldsymbol{P}} \\ {{{\boldsymbol{P}}^{\text{T}}}}&{\bf{0}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\boldsymbol{H}} \\ {\bf{0}} \end{array}} \right] $$ (26) 此时，在点I处计算的插值公式表示为：
$$ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{H}}_I}} \\ {\bf{0}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{A}}_I}}&{{{\boldsymbol{P}}_I}} \\ {{\boldsymbol{P}}_I^{\text{T}}}&{\bf{0}} \end{array}} \right]{\left[ {\begin{array}{*{20}{c}} {\boldsymbol{A}}&{\boldsymbol{P}} \\ {{{\boldsymbol{P}}^{\text{T}}}}&{\bf{0}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\boldsymbol{H}} \\ {\bf{0}} \end{array}} \right] $$ (27) 由此可获取上下限的包络值。
-
为验证上述EMD分解算法以及直观反映其分解步骤和效果，尝试对构造函数进行EMD分解，模拟非平稳信号EMD分解具体流程和结果。构造函数：
$$ {F(n) = 8\cos ( {\dfrac{{2{\text π}n}}{{100}} + \dfrac{{\text π}}{2}} ) + 10{{\text{e}}^{ - \tfrac{n}{{500}}}}\cos ( {\dfrac{{2{\text π}n}}{{70}}} ) + } {\cos ( {\dfrac{{2{\text π}n}}{6}} )} $$ (28) n取值范围为
$0 \leqslant n \leqslant 500$ ，采样点数为500，其对应频率约为5.00 Hz、7.14 Hz、83.00 Hz、此时信号函数图像及其频谱如图3所示。
图 3 构造函数信号波形图
对上述构造函数进行EMD分解，如图4所示。
可以看出，信号模态1分解为高频
$ \cos ( {\dfrac{{2{\text π}n}}{6}} ) $ 分量，而模态2大致为两个较低频率$ 8\cos ( {\dfrac{{2{\text π}n}}{{100}} + \dfrac{{\text π}}{2}} ) 和 10{{\text{e}}^{ - \tfrac{n}{{500}}}}\cos ( {\dfrac{{2{\text π}n}}{{70}}} ) $ 的叠加分量，与常规的带通滤波特性不同，并非完全分解为余弦和余弦阻尼两个模态，案例也凸显出了EMD分解时，同一个模态中包含了多个通带的频率特征，并且关于时域波形具有自适应的特点。
图 4 构造函数EMD分解
模态 3 分量的幅度值较小，频率更低，表现出的特征更接近于信号的趋势性分量，对信号重构影响较小。这样，通过时域自适应分解算法，将信号分解为不同的频率分量，有效分解为不同的特征尺度，便于进一步分析信号。
对该信号进行希尔伯特变换，并解析信号的构造，分解为瞬时频率，绘制时频分布如图5a所示。
图 5 构造函数的EMD分解时频分布图
与图3b对比，可以看出经过EMD分解，成功将拟合信号分解为4个分量，其中，imf1对应原信号83.00 Hz分量；imf2对应原信号7.14 Hz分量；imf3对应原信号5.00 Hz分量；res.为信号残差趋势分量。各分量能够较好地与原信号频域对应，并且获取了各个分量的时频谱，能够进行各个分量的时变频谱分析。
-
为验证BEMD快速算法的合理性及处理效果，建立了二维FDTD模型进行仿真计算，其中，网格区域为
$8\;{\rm{m}} \times 8\;{\rm{m}}$ ，地层介电常数为${\varepsilon _r} = 7$ ，电导率为$\sigma = 1 \times 10^{-4}\;S/m$ ，接收天线和发射天线间距0.2 m，发射频率为80 MHz零阶高斯脉冲，雷达探测深度为8 m，分辨率为0.7 m，目标为半径0.3 m金属圆形目标，时窗选取$T = 2 \times 10^{-7}\;s$ 。通过FDTD正演仿真技术，建立FDTD网格进行正演仿真模拟，对仿真结果进行加噪处理，以分析BEMD对其处理效果。加入具有标准正态分布，即以0为均值、1为标准差的二维高斯噪声，经过分析可得，通过BEMD能够分解B扫数据，更清晰地将目标回波信号反映出来。将imf2、imf3数据进行叠加，可以看出，通过BEMD分解和对其分量的重组，使得目标的回波双曲线更加明显，特征更加突出，提高了B扫数据的信噪比，有助于进一步的雷达信号成像与反演。通过对仿真数据进行分解，验证了BEMD对于井中雷达B扫数据处理的降噪能力。 -
基于理论和仿真计算，对实测数据进行EMD分解，采用的原始数据为油井测试数据，此时并没有对原始数据设定分解信号的基函数，通过信号的自适应包络曲线以及迭代，将信号分解为如图6所示的多个不同特征的imf分量和残差。
其中各IMF分量标注如图，该信号共分解为7个IMF模态分量以及一个单调趋势函数分量res，能看到每一次imf分解完毕，获得的模态分量在时窗范围内是平稳的，其分解得到的特征模态分量从上至下是由高频到低频的顺序进行排列，即imf1是信号的高频主成部分，依次频率减小，res.为频率接近0的残余分量。
模态特征函数分量从上至下是严格由高频到低频的顺序排列，此处并非 imf1频率要高于imf2，而是允许后者局部频率高于前者局部频率，但在同一时间采样点处，频率由高到低排列。如图7所示，框内标注处的分量局部信号，包含有较低频率的局部信号。
图 6 井中雷达信号EMD分解
图 7 EMD分解的局部信号
这也正好反映了EMD分解的局部性强的特性，即相邻的imf可能包含有相同频率范围内的震荡，但在相同的频率范围内的震荡不会出现在同一时窗范围内。
将上述信号进行时频域分解，由于信号的不稳定性，分解之后的时频谱起伏较为剧烈，随时间采样点数增加，频率变化范围较大，频率起伏范围可达300 MHz。通过时频谱可以看出，imf1、imf2、imf3主要包含的是频率为100～300 MHz分量，与井中雷达设计的系统工作频率相符合，因此取 imf1、imf2、imf3这3个分量进行叠加，叠加之后的信号形成B扫图像。与原信号B扫图像对比可得，对原B扫数据按行模态经验分解并叠加imf1、imf2、imf3后，合成的B扫信号能够更多的对地质异常目标进行反映，异常目标更清晰明显，异质构造特征更突出，对分析地层异常提供更合理有效的解释。从分解过程也容易看出，EMD分解还具有线性、完备性的特点。
-
对实测井中雷达数据进行BEMD分解，所选取的数据为中石化榆林大125井约1975～2995 m深井中雷达通道1数据，B扫数据去除背景后通过对其进行BEMD分解，同时叠加imf2、imf3分量，能够很好地剔除原始数据中由于硬件采集模块带来的采样噪点以及其他各部分系统模块带来的系统噪声，减小高频干扰，同时通过分离低频分量，剔除低频干扰。因此，地层波动情况得以更清晰地反映出来，使得回波双曲线更加突出，对于分析地层的不同情况提供更好的信息量支持。
由上述讨论可知，基于EMD分解的井中雷达数据处理通过模态经验分解算法，根据信号本身的时间尺度特性对信号进行分解，有效地克服了傅里叶变换的局限性，更适应于井中雷达信号非线性、非平稳的特点，可以在任何时刻获得信号的频率分布，在时频域上提供更高的频率分辨率。可见，在井中雷达数据分析处理方面基于EMD分解的方法优于小波变换，该方法很好地解决了雷达数据的模式混合问题，保证了分解对噪声的稳定性。针对不同类型储层，特别是缝洞型碳酸盐岩储层这种地层结构复杂的情况，基于EMD分解的方法相对于其他常规信号处理技术，可以获取到更多地质信息。但其模式的分裂问题仍然是一个开放的问题，更重要的是其数学分析并不完整，在理论分析和理解、性能增强、优化等方面仍存在很多问题。
-
本文简要介绍了井中雷达信号模型和数据类型的数学描述以及常规信号预处理方法，对井中雷达信号预处理的理论进行了推导研究。对希尔伯特−黄变换方法进行了详细推导，也进行了基于EMD信号预处理方法的理论研究。通过仿真的方式验证了将EMD的一维方法推广到了实测二维B扫数据方法的正确性、有效性。通过实测数据对各算法进行验证，其中，EMD有效提高了采集原始数据的信噪比；BEMD能够有效剔除高频、低频的干扰，较一维的模态经验分解有较大的改进。结合仿真和实测数据的处理结果可以看出，基于EMD的时频域方法能够有效对井周地层解释提供依据，具备应用到井中雷达数据处理的潜力。
本文研究工作得到了中国石油化工股份有限公司西北油田分公司项目(34400007-20-ZC0607-0050)的支持，在此表示感谢。
EMD-Based Borehole Radar Signal Preprocessing Method
-
摘要: 针对井中雷达数据因环境复杂易受到“污染”而很难提取出有效数据的问题，基于经验模态分解方法(EMD)提出了二维模态经验分解方法(BEMD)，该方法对井中雷达信号预处理的时频域分解算法进行分析，将其拓展至高维进行信号处理。仿真结果表明，二维模态经验分解能够有效剔除高频、低频的干扰，较一维的模态经验分解有较大的改进，验证了该方法具备应用于井中雷达数据处理的潜力。Abstract: Aiming at the problem that the radar data in the well is difficult to extract effective data due to factors such as pollution in the complex environment, this paper proposes a two-dimensional modal empirical decomposition method (BEMD) based on the empirical mode decomposition method (EMD). The time-frequency domain decomposition algorithm of radar signal preprocessing is analyzed, and the two-dimensional radar signals in wells can be processed effectively by extending the signals to high dimensions. The simulation results show that the proposed BEMD can effectively eliminate the interference of high frequency and low frequency and achieves a great improvement compared with the one-dimensional modal empirical decomposition, which verifies that the method has the potential to be applied to the processing of radar data in Wells..
-
Key words:
- borehole radar /
- EMD algorithm /
- finite difference time domain /
- signal preprocessing
-
[1] PERSICO R. Introduction to ground penetrating radar: Inverse scattering and data processing[M]. [S.l.]: John Wiley & Sons, 2014. [2] 陈建胜, 陈从新. 钻孔雷达技术的发展和现状[J]. 地球物理学进展, 2008, 23(5): 1634-1640. CHEN J S, CHEN C X. Development and current situation of borehole radar technology[J]. Advances in Geophysics, 2008, 23(5): 1634-1640. [3] NEAL A. Ground-penetrating radar and its use in sedimentology: Principles, problems and pro-gress[J]. Earth Science Reviews, 2004, 66(3-4): 261-330. doi: 10.1016/j.earscirev.2004.01.004 [4] OLSSON O, FALK L, FORSLUND O, et al. Borehole radar applied to the characterization of hydrau-lically conductive fracture zones in crystalline rock[J]. Geophysical Prospecting, 1992, 40(2): 109-142. doi: 10.1111/j.1365-2478.1992.tb00367.x [5] RANNEY K, STANTON B, NGUYEN L, et al. Borehole radar performance characteristics and applications for underground change detection[C]//2006 IEEE Conference on Radar. [S.l.]: IEEE, 2006: 7. [6] HOLSER W T, BROWN R J S, ROBERTS F A, et al. Radar logging of a salt dome[J]. Geophysics, 1972, 37(5): 889-906. doi: 10.1190/1.1440307 [7] 马春光. 瞬态脉冲雷达成像测井及实验研究[D]. 成都: 电子科技大学, 2015. MA C G. Transient pulse radar imaging logging and experimental research[D]. Chengdu: University of Electronic Science and Technology of China, 2015. [8] TRICKETT J, MASON I, STEVENSON F. Borehole radar at an underground Ventersdorp reef site[C]//The 6th South African Geophysical Association/SEG Biennial Conference. Cape Town: [s.n.], 1999: cp-221-00027. [9] TIRÉN S W C. Borehole radar measurements aid structure geological interpretations[J]. Journal of Applied Geophysics, 2000, 43(2): 227-237. [10] HALABE U B, CHEN H, BHANDARKAR V, et al. Detection of sub-surface anomalies in concretebridge decks using ground penetrating radar[J]. ACI Materials Journal, 1997, 94(5): 396-408. [11] OSWALD B, BENEDICKTER H, BACHTOLD W, et al. A single-rod probe for time domain reflectometry measurements of the water content[J]. Vadose Zone Journal, 2004, 3(4): 1152-1159. doi: 10.2136/vzj2004.1152 [12] GRASMUECK M. 3-D ground-penetrating radar applied to fracture imaging in gneiss[J]. Geophysics, 2012, 61(4): 1050. [13] EL-MAHALLAWY M S, HASHIM M. Material classification of underground utilities from GPR images using DCT-based SVM approach[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 10(6): 1542-1546. doi: 10.1109/LGRS.2013.2261796 [14] ALEJOS A V. Understanding the design of anti-dispersive filtering for propagation of UWB microwave signals in dispersive soils[J]. IEEE Geoscience and Remote Sensing Letters, 2013, 11(1): 14-18. [15] 张先武, 高云泽, 方广有. 消除探地雷达数据的子波衰减和频散的反滤波方法[J]. 地球物理学报, 2014, 57(3): 932-938. doi: 10.6038/cjg20140322 ZHANG X W, GAO Y Z, FANG G Y. Inverse filtering method to eliminate wavelet attenuation and dispersion of ground penetrating radar data[J]. Acta Geophysics, 2014, 57(3): 932-938. doi: 10.6038/cjg20140322 [16] 王立凯. 希尔伯特−黄变换在雷达测井中的应用研究[D]. 成都: 电子科技大学, 2016. WANG L K. Application of Hilbert-Huang transform in radar logging[D]. Chengdu: University of Electronic Science and Technology of China, 2016. [17] WU Z, HUANG N E. Ensemble empirical mode decomposition: A noise-assisted data analysis method[J]. Advances in Adaptive Data Analysis, 2009, 1(1): 1-41. doi: 10.1142/S1793536909000047 [18] WU Z, HUANG N E, CHEN X. The multi-dimensional ensemble empirical mode decomposition method[J]. Advances in Adaptive Data Analysis, 2009, 1(3): 339-372. doi: 10.1142/S1793536909000187 [19] GAN L, ZHOU L, YOU X, et al. The instantaneous frequency extraction of GPR B-scan data based on HHT method[C]//2012 International Conference on Machine Learning and Cybernetics. [S.l.]: IEEE, 2012: 982-985. [20] LI J, LIU L, ZENG Z, et al. Simulation and signal processing of UWB radar for human detection in complex environment[C]//2012 the 14th International Conference on Ground Penetrating Radar (GPR). [S.l.]: IEEE, 2012: 209-213. [21] ADDISON A D, BATTISTA B M, KNAPP C. Improved hydrogeophysical parameter estimation from empirical mode decomposition processed ground penetrating radar data[J]. Journal of Environmental & Engineering Geophysics, 2009, 14(4): 171-178. [22] CHEN C S, JENG Y. Nonlinear data processing method for the signal enhancement of GPR data[J]. Journal of Applied Geophysics, 2011, 75(1): 113-123. doi: 10.1016/j.jappgeo.2011.06.017 [23] JENG Y, CHEN C S. A nonlinear method of removing harmonic noise in geophysical data[J]. Nonlinear Processes in Geophysics, 2011, 18(3): 367. doi: 10.5194/npg-18-367-2011 [24] WU X, SENALIK C A, WACKER J, et al. Object detection of ground-penetrating radar signals using empirical mode decomposition and dynamic time warping[J]. Forests, 2020, 11(2): 230. doi: 10.3390/f11020230 [25] CHEN Y C, KAO S P, HSU Y C, et al. Automatic measurement and computation of stream cross-sectional area using ground penetrating radar with an adaptive filter and empirical mode decomposition[J]. Journal of Hydrology, 2021, 600: 126665. doi: 10.1016/j.jhydrol.2021.126665 -
点击查看大图
图(7)
计量
- 文章访问数: 3496
- HTML全文浏览量: 1008
- PDF下载量: 61
- 被引次数: 0
