1 引言

在过去的几十年间，探地雷达在理论研究、技术、设备制造工艺、以及应用方面都取得了重大进展.作为一种非常有效的地球物理方法技术，探地雷达被广泛应用在地下水污染探测评价、地质工程、沉积学、冰川学和考古学中(Bristow and Jol, 2003; Daniels, 2004; Jol, 2009; 傅磊等, 2014).

探地雷达系统包括发射天线和接收天线，以及控制信号收发和数据存储的控制系统.探地雷达使用宽带天线，不同天线具有不同的频率范围，同时也具有不同的探测深度和分辨率.探地雷达采用较低频率进行探测时，可获得的探测深度较大，但分辨率相对较低；反之采用较高频率的电磁波进行探测时，探测深度浅，但能获得相对较高的分辨率.

钻孔雷达是探地雷达中的一种，由于天线放置在钻孔中，它可以探测和评估地下深部区域地质结构，并且具有很高的分辨率.钻孔雷达的应用包括地下结构探测、特征描述、地质研究，以及地下设施探测等.目前的钻孔雷达使用的天线频率一般低于100 MHz，在结晶岩中的穿透深度为20 m或更远(Sandberg et al., 1991; Olsson et al., 1992；Sato and Takayama, 2007).

传统的钻孔雷达发射天线和接收天线都是全方位的，因此，只能探测目标地质体深度，以及其对井孔的距离，并不能对其方位进行探测.目前，对地下介质进行三维立体成像仍然是一件非常困难的事情.为了解决这个问题，发展了定向钻孔雷达技术(刘四新等，2005).

定向钻孔雷达不仅可以探测目标的深度和距井孔距离，而且可以探测其方位.定向钻孔雷达的实现大致有两种方式：其一使用阵列天线接收的方式进行测量，主要使用均匀圆阵，根据不同的波至时间提取方位信息；其二采用定向发射天线的方式，天线向井周定向发射电磁波，目标体的方位位于天线的发射方向.在常规的钻孔雷达系统中，发射天线向四周同等程度地发射电磁波，并没有特别的方向性.但是，在定向发射天线的钻孔雷达系统中，需要使得发射天线具有明确的方向性，这样才能通过反射波准确地判定目标地质体的方位(曾昭发等，2006).

美国加利福尼亚大学劳伦兹利物莫实验室的Lytle和Laine(1978)，最早开发定向钻孔雷达的发射天线.这种雷达的天线包括一个点辐射元，其位于高介电常数介质的偏心位置.由于发射天线偏心所造成的影响，天线在转动的过程中，向各个方向所发射的电磁波并不是完全相同的，而是有一定方向性的.

荷兰的相关部门研制了一种定向钻孔雷达仪器，其采用旋转的定向发射天线.这种雷达的天线系统包括一个普通的发射天线和一个屏蔽系统，发射天线向各个方向发射强度相等的电磁波，屏蔽系统是用特殊的材料制作的，可以用来吸收发射天线的电磁辐射.因此，特定方向上的电磁波被屏蔽掉，这时的发射天线便获得了方向性.当发射电磁波具有较高的频率时，它所判定的方位是非常准确的，但是它的探测范围是比较小的.因此，在研制的过程中，既要保证定向钻孔雷达判定方位的准确性，又要保证该雷达具有一定的探测范围.

定向接收型的钻孔雷达较为常见.作为定向接收天线的一种，所使用的天线是两个正交的环状天线.雷达所发射的电磁波经过目标地质体反射回来，并且入射到接收天线.当入射波和天线所在的平面处于垂直关系时，天线接收不到任何反射波信息.当入射波和天线所在的平面处于同一个平面上时，天线接收到的反射波的信号幅度是最大的.但是，这种定向钻孔雷达技术也存在着明显的缺点，环状天线的大小受到了井眼大小的限制，影响了信号检测.另外，只能将目标体定位为两个相距180°的方向，其具体的方位的确定还需其他辅助信息.

由多个小天线组成的圆环天线阵列成为定向钻孔雷达定向接收技术的研究热点(Ebihara and Yamamoto, 2006; Ebihara and Inoue, 2009; Ebihara et al., 2011, 2013, 2015).Ebihara等(2000)将MUSIC(Multiple Signal Classification，多信号分类)算法引入定向钻孔雷达系统，并取得了较好的效果，同时确定了MUSIC算法在低信噪比情况下的局限性.Takayama和Sato(2007)将反正切法引入定向钻孔雷达系统，Ebihara等(2013)提出了使用计算残差来确定目标方位的算法.可以看出，随着定向钻孔雷达技术的进一步发展，定向钻孔雷达系统和方位估计算法均有了较大的发展, 但方位确定和三维成像算法依然需要进一步解决和完善.

本文研究的定向钻孔雷达为一个发射天线和四个接收天线(如图 1)，四天线分布于一个垂直于井眼的圆环位置上，呈等角度分布.为了分析方便，我们把四个天线分别指定为东、南、西、北天线，顺时针旋转.当北天线方向和地理正北一致时，其余天线分别和东、西、南一致，否则会有偏差角，这可以通过测量北天线相对于地理正北的夹角来校正.所以，计算目标的方位只需要计算相对于四天线的方位即可.我们将四个接收天线所接收到的波形进行对比分析，利用相应的方位估计算法，确定目标地质体的具体方位角，进而利用方位角对地质体进行三维成像.

本文提出并实现了基于BS-MUSIC(Beam-Space based MUSIC的简称)的钻孔雷达方位确定方法，并实现了井旁目标的三维成像方法.BS-MUSIC算法是在MUSIC算法基础上的改进，将四个接收天线所接收到的信号合成两个波束，达到提升运算速度的目的.利用数值模拟结果进行验证，取得了非常好的效果.

2 BS-MUSIC算法

BS-MUSIC算法和MUSIC算法都是建立在均匀圆阵理论上的(王永良等，2004).

2.1 均匀圆阵的一般理论

首先探讨入射波和阵列的关系.均匀圆阵如图 2所示.若天线阵由M个全向天线组成，有K个信源，θ, φ分别为球坐标系中目标的仰角和方位角，第i个信号为：

(1)

式中，z_i(t)为第i个信号的复包络，包含信号信息.e^jω₀t为空间信号的载波，ω₀为角频率.则整个天线阵接收到的信号为：

(2)

式中，A=[a₀ … a_K－1]为阵列流形，a_i=[e^{－jω₀τ_1t} … e^{－jω₀τ_Mt}]^T为信号i的方向向量，τ为相对于圆心的时间延迟，S(t)=[s₀(t) … s_M－1(t)]^T为信号矩阵，N(t)=[n₁(t) … n_M(t)]^T为加性噪声矩阵，[·]^T表示矩阵转置.

均匀圆阵相对于波达方向(Direction Of Arrival, DOA)为φ的信号的方向向量为：

(3)

式中γ_m=2πm/M，ξ=k₀Rsinθ，k₀=2π/λ₀，λ₀为波长，R为圆阵半径.

2.2 基于BS-MUSIC算法理论

MUSIC算法利用了信号子空间和噪声子空间的正交性，构造空间谱函数，通过谱峰搜索，检测信号的波达方向，峰值对应目标方位.

MUSIC算法是建立在以下假设基础上的：

(1) 阵元间距不小于处理最高频率信号波长的二分之一；

(2) 处理器的噪声为加性高斯分布；

(3) 信号源少于阵元数，信号取样数大于阵元数；

(4) 信源必须为均匀分布的随机数.

BS-MUSIC是MUSIC算法的改进，主要是先将空间阵元通过变换合成一个或几个波束，再利用合成的波束数据，进行DOA估计，其原理如图 3所示.相对于传统的MUSIC算法，基于BS-MUSIC算法有以下明显好处：一是降低计算量；二是提高空间谱估计方法的稳定性；三是降低系统的复杂性.

可以看出，波束空间处理需要通过阵元合成一定数量的波束通道作为数据的接收通道，这和阵元空间处理中每个阵元对应一个接收通道完全不相同.如果我们将基于特征分解的信号子空间算法的运算量设为O(M³)，对于大阵列、小信号源数的场合，这类算法很难实时实现.显然，波束空间方法在这种场合下，能有效降低算法的计算量，如果合成波束的维数L，且K＜L＜M(其中K为信源的个数，M为阵元的个数)，则DOA估计的计算量由O(M³)降为O(L³)，所以，在合成的波束数L远小于阵元数M时算法的运算量显著下降.同样，就系统设计而言，一个数据接收通道就对应一个接收机，显然，数据接收通道的减少可以大大减少系统复杂度.

对于半径为λ/2的均匀圆阵，θ=π/2时，阵列的方向向量

(4)

可定义如下波束形成矩阵：

(5)

可以看出，式(5)的每一列表示波束主瓣指向u=sin(2kπ/M)的波束形成器，其中k=0, 1, …, M-1，则式(5)共有M个波束形成器，各相邻主瓣指向的间隔为Δ=2π/M.

由此，我们可以构造波束形成矩阵，并进行归一化加权

(6)

由于波束输出由M×1向量变为L×1向量，对应的协方差矩阵为

(7)

其中B=T^TA=[T^Ta(φ₁) … T^Ta(φ_N)]=[b(φ₁) … b(φ_N)]，R_S=XX^T.对式(7)进行特征分解得到特征值，其中包括一个大的特征值和若干较小特征值，去除最大特征值后较小特征值对应的特征向量组成噪声子空间，可得到BS-MUSIC算法的空间谱.

(8)

根据空间谱的峰值位置可以判定信号的方位.

对于BS-MUSIC算法，需要说明的是：

(1) BS-MUSIC算法的实质是先将阵列形成几个波束，换言之，就是将阵列接收的数据，变换到波束空间的数据，再利用波束空间的数据，进行DOA估计.

(2) BS-MUSIC算法在利用MUSIC算法进行搜索时，注意T的维数必须等于形成的波束数，而且，在对波束空间的数据进行分解时，矩阵的维数也由原来的阵元数降低为波束数.

(3) 波束形成矩阵T的选择决定波束空间算法的性能.

3 基于BS-MUSIC的定向钻孔雷达的成像算法

使用BS-MUSIC算法和MUSIC算法对定向钻孔雷达数据进行处理前，需进行以下二步处理，以满足算法要求.

(1) 复信号转换.MUSIC算法是针对于时间域的复信号进行处理，而定向钻孔雷达所接收到的数据为实信号，所以需要将实信号转化为复信号.观察信号组成可以发现，z_i(t)e^jω₀τ，其中z_i(t)是信号部分，通过随机发生器产生.由于窄带假设，四个接收器所接收到的信号部分是一样的，只是存在不同的延迟，因此一切确定方位算法的核心都是对于波形到达的先后次序进行处理，以此来估计目标的精确方位，我们将四个天线分别置于东南西北四个位置，正是为了得到波形到达的先后次序，其中j为虚部，ω₀由实际信号的主频获得，而τ为相对于圆阵中心的时间延迟.

如图 4所示，圆阵中心为O，东、西、南、北四个接收天线分别位于E、W、N、S位置.入射波穿过圆阵中心，可以明显看出OA=OB，也就是东向接收天线所接收到的信号相对于圆阵中心的时间延迟大小和西向接收天线是一样的，即τ_E=-τ_W，这样，只要求出东向接收天线相对于西向接收天线所接收到信号的时间延迟，除以2，便可以得到τ_E、τ_W的具体数值.我们采用残差法计算时间延迟.同理，对于南向接收天线和北向接收天线也是适用的.这样，四个方向接收天线的时间延迟就确定下来.由于空间谱方法仅对均值为0的随机信号适用，因此，通过随机函数，产生一组随机信号，将时间延迟与随机信号结合，实信号就成功转化为复随机信号.

(2) 半波长圆阵信号计算.对于BS-MUSIC算法而言，当圆阵半径为λ/2时，算法的效果比较好.实际中由于钻孔半径很小，不满足这一要求，需对信号进行转换，计算出半波长圆阵等效的信号.如图 5所示，其中内圈部分的四个端点，即E、N、S、W为四个实际的接收天线位置，外圈部分为假想的半径为λ/2圆阵，四个端点TE、TN、TS、TW为假想的四个接收天线，假想天线相对于圆心的时间延迟可以通过四个实际接收天线的时间延迟求出来.

完成以上校正之后，我们就可以用MUSIC算法和基于BS-MUSIC算法进行处理.在BS-MUSIC算法中，采用两个合成波束进行运算.逐道进行处理.在处理每一道数据时候，先选择一个时间窗口，确保在时间窗口内部只有一个波形占据主导地位，并且当信号强度大于阈值才进行处理，否则被当成噪声忽略掉.另外，只有波形到达先后可以被区分的情况下才进行处理，如果四个波形同时到达，则不作判断，将这块信号忽略掉.计算时，我们先使用之前得到的时间延迟和一个均匀分布的随机信号进行相乘，得到新的数据，然后使用MUSIC算法和基于BS-MUSIC算法进行计算.

基于BS-MUSIC算法的流程总结如下：

(a) 输入计算所需参数：天线收发距、天线主频、圆阵半径、电磁波在介质中的传播速度、信源个数、时间采样间隔等.

(b) 对信号进行必要的预处理，包括去除直达波和去噪等，以消除干扰.

(c) 利用经过预处理的数据进行方位确定.首先计算东南西北四个接收天线相对于圆阵中心的时间延迟，然后换算成半径为半波长的均匀圆阵相对于圆阵中心的时间延迟.将所求的时间延迟和随机信号结合得到X矩阵，求取波束变换矩阵B，并将波束变换矩阵与X相乘得到协方差矩阵R_yy，将协方差求取特征值和特征向量，再利用BS-MUSIC算法计算得到方位角.由于采用逐点处理，得到的方位信息也是一个时间序列.

有了方位信息后，可以利用得到的数据进行三维成像，三维成像过程总结如下：

(a) 对于每个深度点，输入方位角和四个接收天线接收到的信号，并对数据进行增益处理.

(b) 四个接收天线所接收到的信号的平均信号作为该深度点的信号，利用计算得到的方位角数据把信号分配到各个方位，这样在各个方位都有一个时间序列信号.方位间隔呈等角度分布，比如，10°或5°间隔.这样在每个深度点，得到一个横向切片.对各深度点都做同样处理，形成一个三维数据体，进而得到各深度对应的横向切片和各方位对应的纵向切片，分别对应10°，20°，…，360°的信号(以10°间隔为例).

(c) 对各个方位角度对应的纵向切片数据进行偏移处理，这时三维数据即为三维图像数据.这里采用绕射叠加的偏移算法.为了便于显示，数据需从柱坐标系转换为直角坐标系.

4 合成数据算例

我们利用GprMax(Giannopoulos，2005)软件合成定向钻孔雷达数据，检验提出的算法.模型包含一个裂缝和一个球形空洞，裂缝在井的正东方向，球形空洞在井西南方向，如图 6所示.模拟中没有考虑井眼的影响.裂缝深度4 m到8 m，宽3 m，平行于井，中心距井2.5 m；球深度为10 m，球心位于井西南方向7 m，球半径为1 m.围岩采用花岗岩，裂缝和球形空洞内填充物为水.花岗岩物理参数选取如下：相对介电常数为5，电导率为0.0007，相对磁导率为1.水的参数选取如下：相对介电常数为80，电导率为0.7，相对磁导率为1.模型被离散为均匀网格，网格大小dx=dy=dz=0.1 m.时窗大小为320 ns.

本次模拟采用一个发射天线和四个接收天线，四个接收天线分别位于东南西北四个方向.发射天线位于井眼中间，四个接收天线到井孔中心的距离为0.1 m，收发距为2.5 m，天线主频为100 MHz，发射天线和接收天线同时从井底向上移动，移动步长为0.1 m，移动次数为100次.相应的合成数据波形如图 7所示，这里采用每十道显示一道的方式，以便能看清波形.从波形数据可明显看出直达波，和来自两目标的反射波，且各天线的反射波至时间稍有差别，这些差别为定向探测提供了可能.

对于钻孔雷达所采集的数据，首先要去除直达波等干扰成分.接下来进行复信号转换，再计算半波长圆阵列对应信号.完成转化之后，采用MUSIC算法和BS-MUSIC算法进行处理，处理结果使用滤波的方法滤掉由于计算误差产生的毛刺.

处理过程采用64位电脑处理器：Intel(R) Core(TM) i5-3230M CPU@2.60 GHz, 内存8.0 GB.对于同样数据，经典MUSIC算法耗费14.638439 min，而BS-MUSIC算法仅仅耗费6.037945 min，可以看出，速度有了非常明显的提升.

选取第75道数据处理的中间结果进行展示.图 8a和8b分别为MUSIC算法和基于BS-MUSIC算法计算的空间谱，图的横坐标代表的是走时，纵坐标代表的是每个方位角对应的空间谱，每个走时均对应360个方位角，每个方位角都有一个空间谱数值，其中空间谱数值最大的方位角，即为该走时的波达方向，并且最大值与其他方位角对应空间谱的数值差的绝对值越大，说明该方法的效果越理想，在图中出现大量空间谱为0的区域，这是因为在对波形进行计算的过程中，该点的振幅低于阈值，则当成噪声，不进行方位识别，并赋值为0.从图 8a中可以明显看出，在90°和225°处有两个最大值，对应目标方位，和模型符合，从图 7b也可以明显看出90°和225°两个最大值，其中MUSIC算法的目标方位角对应的空间谱数值和其他方位角对应的空间谱数值差距并不大，说明MUSIC算法的性能相对较差，而BS-MUSIC算法的区分度更加明显，再次说明了BS-MUSIC算法的性能相对较好.这两种算法的识别方位的效果会随着信噪比的降低而逐渐变差，但是这个过程中BS-MUSIC算法的效果会优于MUSIC算法，因为从目前的情况来看，BS-MUSIC算法的区分度更高一些，因此在取最大值的时候效果会更好一些.

图 8c和8d分别为MUSIC算法和基于BS-MUSIC算法计算的方位，可以看出，这两种算法的处理精度都是非常高，误差并没有超过1°的情况.

结合以上计算得到方位角和合成信号，可以计算三维数据体，进而形成横向切片和纵向切片数据.两种算法得到的模型的横向切片如图 8e和8f所示，从图中可以清楚看出正东方向和西南方向存在一个目标地质体，这与实际情况是符合的.两种算法得到的模型的纵向切片如图 9和10所示，从这两图中可以明显看出，90°方向有一个信号，西南225°方向也看到明确信号，这也和模型是符合的.从横向和纵向切面来看，方位确定取得很好效果.

使用经典MUSIC算法和BS-MUSIC算法得到的三维成像结果如图 11a和11b所示.可以看出，偏移处理使绕射波得到归位，取得非常好的成像效果.从图 11可以明确看出，正东方向存在一个裂缝，由于定向钻孔雷达仅仅在垂直方向上移动，因此，可以在垂直方向上有效判定目标地质体的规模，在水平方向上由于缺少数据，因此在水平方向上无法确定目标规模形状等，但是可以根据走时计算出目标距井的距离.

5 结论

研究的定向钻孔雷达包含一个发射天线和四个接收天线，通过所接收到波形到达的先后次序，可以有效地确定目标地质体的精确方位.本文提出并实现了基于BS-MUSIC算法的目标方位确定方法，在此基础上实现了钻孔雷达三维成像方法.使用GprMax合成数据，对基于MUSIC算法和基于BS-MUSIC算法的结果进行比较，可以看出两者均可以有效判断目标地质体的精确方位，并且在三维成像中两者均取得了非常好的效果.同时，我们发现，BS-MUSIC算法在运算速度上明显优于MUSIC算法，非常适合定向钻孔三维成像.