当地震波通过应力变化区时其传播速度也要发生变化，由于P波和S波对地壳介质物理参数的响应不同，因此P波和S波速度的变化不相同，两者的比值就发生变化(国家地震局科技监测司，1995；彭朝勇等，2013a；Peng et al., 2014a，b).波速是地震波的一种运动学特征，较大地震发生前，孕育区内岩石的波速可能出现某些变化(即波速异常)，如何识别地震波波速异常并研究其特点，对于探索地震预报途径及震源孕育过程有着重要的实际意义(陈章立等，2007).新丰江6.1级水库地震和我国成功预报的海城地震在地震前后都存在波速变化(冯锐等，1976；冯锐，1977).用多台法和固定台站方法对云南地区的永胜、大姚、姚安、施甸地震序列的波速比计算发现主震发生后，在下一次强余震发生前，波速比出现了趋势性下降—回升—发震的过程(黎明晓和刘杰，2006).我国大陆灾害性的地震几乎全都发生在震源深度为5～25 km的地壳层中，因此对地壳层物理参数变化的动态监测是实现地震预报的一个行之有效的办法(庄灿涛，2002).由于天然地震在时间、空间和强度上的不确定，因此无法通过天然地震来对地下结构进行连续观测，以了解地下介质的波速变化动态.人工爆破在时间、空间和强度上可以控制，但由于花费较高、破坏性较大和不具有重复性，也不适合用作对地下介质波速变化的连续动态监测.气枪震源用水作为传播介质，它对水源的深度和体积要求高，在高寒地区水结冰后就不能工作了.

利用绿色、环保的人工震源，主动向地下发送地震波构建地震波雷达，实现大范围的地下探测是地球科学的前沿课题(陈颙等，2006).地震波雷达(也叫精密可控震源)是20世纪末发展起来的一种主动探测新技术，通过精密控制振动系统，对目标区域较长时间发射低能量的信号，用高灵敏度低噪声的地震数据接收设备和现代地震数据处理方法，可以把这种长时间低能量的记录数据转变成可以用于反演地下结构及物性参数的形式.地震波雷达运行时间与GPS同步，精确到1 ms，运行地点由GPS精确测定，占地面积为9 m²，相对天然地震而言，它是一个很准确的“点”源.地震波雷达是机电设备，控制精确、功耗低、噪音小、重复性高、无破坏性、可移动使用，是一种绿色环保的震源.俄罗斯在西伯利亚，距地震波雷达300 km外提取到了信号(Alekssev et al., 2005)，日本用单频地震波雷达得到地下波速变化(Ikuta and Yamaoka， 2004).我们国家应用地震波雷达才刚刚起步，只有有限的几位学者开展了该方面的研究，如王洪体等(2009)用输出力为40 t的地震波雷达在河北沽源计算出了地震波的走时变化，李志伟等(Li et al., 2013；刘希康等，2013)在广东新丰江提取到180 km以外输出力为40 t的地震波雷达发出信号的走时，杨微等(2010)用输出力为10 t的地震波雷达(移动式)测量数据监测到了四川绵竹M5.6地震前后波速变化.高精度组合式轻便小型可控震源弥补了炸药震源在地表条件差以及城市环保区无法施工的缺陷，在中浅层勘探中取得了较好的效果(赵春蕾等，2013).常用的提取信号延迟方法有各自的优点和不足，互相关方法处理结果稳定性好．能较好地突出高能量频率成分，但不适用于在频谱特征里能量相差较大的窄带信号分析；精密可控震源的信号与环境背景噪声在频谱上有重叠，短时相关法会放大噪声的作用，其效果较相干和反褶积法显得略差一些；信噪比对相干和反褶积法处理结果的稳定性影响较大，相干和反褶积法适合于具有一定信噪比的精密可控震源信号检测分析(杨微等，2013).

地震波雷达主要由两台伺服电机、两个机械特性完全一样的偏心轮和一套控制电路组成(如图 1所示).两台伺服电机分别带动偏心轮相向转动而向地面输出一个垂直向下的合力，通过精确控制两台伺服电机的转动频率而实现连续向地下发射线性调频信号，输出力的表达式(王洪体等，2009):

图 1

Fig. 1

图 1 地震波雷达图 (a)地震波雷达外观结构；(b)地震波雷达运行示意图.Fig. 1 Seismic wave radar (a)Appearance of the SWR structure;(b)SWR running diagram.

假设线性调频的起始频率为f_start，结束频率为f_end，调频周期为T，单位为s，则调频斜率为f_k：

通过编码器和控制电路来控制地震波雷达的调频周期T、运行时间以及运行的开始频率f_start和结束频率f_end，在指定地点和时间发射出指定调频斜率的线性调频信号.线性调频信号是一种典型的非平稳信号，它的频率随时间变化而线性变化，具有很好的时间-频率聚集性，被广泛用于雷达、声纳、通信、地震勘探等领域(张贤达和保铮，1998).

地震波雷达从每小时的第300 s开始启动，在60 s内，它的工作频率由0 Hz线性调频上升至4 Hz. 从360～660 s，它的工作频率从4 Hz线性调频上升至10 Hz，调频时间为300 s，调频范围为6 Hz，根据 公式(2)可知线性调频斜率f_k为0.02 Hz/s.从660～960 s，它的工作频率从10 Hz线性调频下降至4 Hz，调频时间同样也为300 s，调频范围为6 Hz，根据公式(2)可知线性调频斜率f_k为-0.02 Hz/s，完成一个线性调频的上升和下降周期.如此重复5个周期后，从第3360～3420 s工作频率从4 Hz线性调频下降至0 Hz，停止运行(如图 2所示).

图 2

Fig. 2

图 2 地震波雷达的调频模式Fig. 2 Frequency modulation model of SWR

项目期间，在河北沽源县安装了1台地震波雷达，为了监测它发出的线性调频信号，在距地震波雷达10 m处(称之为近场监测点)安装了一套FSS-3M三分向短周期地震计和EDAS-24IP三通道地震数据采集记录器(Peng et al., 2013b；彭朝勇等，2015)，由GPS给地震数据采集记录器授时，授时精度1 ms.为了检测地震波雷达发射出的经地下介质传播后的信号，在沿地震波雷达西南方向布设一条由30套仪器组成的测线(称之为远场检测点)，第一个检测点距地震波雷达约5 km，每个检测点之间相 隔约5 km(如图 3所示).同样，每一个检测点都安装了一套FSS-3M地震计和EDAS-24IP地震数据采集记录器.

图 3

Fig. 3

图 3 监测点和检测点分布图(王洪体等，2009) 五角星为地震波雷达和监测点，同心圆为检测点，虚线为H-21第二炮测线.Fig. 3 Distribution map of detecting and monitoring sites(Wang et al., 2009) Asterisk: SWR and monitoring sites，concentric circles: detecting sites，the dashed line: profile of shot 2 from H-21.

地震波雷达每天分两个时间段运行：9 ∶ 00—19 ∶ 00和21 ∶ 00—7 ∶ 00，从2008年5月9日至5月18日的10天时间里，共接收了126个小时的数据，地震波雷达信号近场监测点和远场检测点的数据采样率都是200 sps.

为了检验地震波雷达信号是否在指定的时间发射指定频率的线性调频信号，对近场监测数据做时间-频率分析.时间-频率分布的主要任务是描述信号的频谱含量是怎样随时间变化的，目的是为了建立一种分布，能在时间和频率上同时表示信号的能量或者强度，得到这种分布后，就可以对各种信号进行分析、处理，提取信号所包含的特征信息，或者综合得到具有期望的时频特征分布的信号(胡昌华等，2002).由于地震波雷达发射信号的频率和输出力都随时间的变化而线性变化，是一种的非平稳信号线性调频信号，因此非常适合用时间—频率联合分析.

对于给定的非平稳信号s(t)，短时傅里叶变换定义为(唐向宏等，2008)

地震波雷达在一个线性调频上升期连续300 s向地下输出频率从4 Hz到10 Hz线性调频信号，由于信号能量衰减和地下介质吸收，只有一部分频率的信号到达远场检测点.我们的目的是要知道这些信号到达远场检测点的走时，因此需要把这信号的能量累积起来.为了累计这300 s信号的全部能量E，先把能量聚集在线性调频直线上.从最佳展现线性调频信号的频率调制规律这一意义上讲，Wigner-Ville分布(WVD)具有理想的时频聚集性(张贤达和保铮，1998).WVD是分析非平稳信号时变的重要工具，它的重要意义之一就是具有明确的物理意义，可以看作信号能量在时域和频域中的分布(胡昌华等，2002).信号s(t)的WVD分布可以表示成(张贤达等，1998)：

线性调频信号的WVD为沿线性调频直线分布的冲击线谱，即时间-频率分布的幅度值集中出现在表示信号瞬时频率变化的线性调频直线上.地震波雷达连续向地面发射线性调频信号，近场监测点记录信号在WVD上呈现一条线性调频直线(如图 4所示)，不同颜色亮度表示不同的信号能量，颜色亮度越高，能量越强.远场检测点信号由于经过地下介质不同路径的传播后，有直达波、反射波和折射波等的出现，因此在WVD上呈现多条具有相同倾角的线性调频直线，不同的直线上聚集了相应地震波能量，这是地震波雷达信号检测的基础.

图 4

Fig. 4

图 4 监测点信号的WVDFig. 4 WVD of the monitoring site signal

为了检测地震波雷达发出信号到达远场检测点的时间，把信号能量沿线性调频直线累积起来.用Hough变换，在时间-频率分布中对线性调频信号的检测和估计可以转换成检测和估计图像中的直线.对于检测图像中的直线，采用极坐标表示形式比采用直角坐标表示形式更合适(胡昌华等，2002)，直线的极坐标形式：

把近场监测点记录到的线性调频信号经过WVD变换后，再进行Hough变换，在三维极坐标(ρ，θ，E)上出现一个波峰(如图 5所示)，波峰值就是聚集在该条直线上的全部能量E.波峰的坐标(ρ，θ)提供了线性调频信号参数估计值，从极坐标变换到直角坐标，就可以估计出地震波雷达调频的斜率和开始调频频率.

图 5

Fig. 5

图 5 监测点信号的Hough变换Fig. 5 Hough transfer of the monitoring site signal

在远场检测点记录到的信号，由于经过地下介质不同路径的传播后，出现了直达波、反射波和折射波等波形，通过Wigner-Hough变换，把聚集了不同地震波能量的线性调频直线在三维极坐标内累积成了不同的地震波峰，波峰对应的时间就是地震波到达的相对时间.

地震波雷达从4 Hz运行300 s后工作频率为10 Hz，这期间它的调频率为0.02 Hz/s，可以计算出线性调频直线的倾角θ.对信号进行WVD变换后，信号能量都聚集在倾角为θ的线性调频直线上，再进行Hough变换后，在三维极坐标内角度为θ的信号提取出来，把近场监测点波峰对应的时间作为地震波雷达发射信号开始的参考时间，计算出每个远场检测点的线性调频信号到达时间，这就是通过能量累计法来实现对地震波雷达信号的检测方法.

近场监测点信号的信噪比高(如图 6所示)，从时间-振幅图中可以看到有五大组波形的时间和振幅变化趋势与地震波雷达的调频时间一致.用短时傅里叶变换对它进行时间-频率分析，从时间-频率 图中可以看到5组从4 Hz到10 Hz再到4 Hz的线性调频直线，信号的时间和频率与地震波雷达控制的运行时间和频率完全一致，同时也产生了多次谐波.

图 6

Fig. 6

图 6 监测点信号图 (a)时间-振幅；(b)时间-频率.Fig. 6 Signals of the monitoring sites (a)Time-Amplitude，(b)Time-Frequency.

为了监测地下介质波速的变化情况，要用地震波雷达长期向地下发射信号，地震波雷达运行的重复性显得尤为重要.信号相关为波形相似的度量，通常用相关系数来表示，相关系数越大，则表明波形就越相似，如果两个信号的相关系数为1，则表明这两个波形完全相同.在本次实验中，地震波雷达在每个小时的固定时间开始启动，运行频率为4～10～4 Hz，固定运行5个上升周期和5个下降周期(如图 2所示)，因此用4～10 Hz的带通滤波器对全部监测点数据进行滤波后选择1 h的数据作为参考信号，再与其他小时数据进行互相关分析.相关的结果显示在本次实验的126 h运行结果中，除了3 h有机械故障外，其余的123 h的数据的相关系数大于0.999(如图 7所示)，表明地震波雷达发射信号的相似度非常高的.在长期的运行中，为了避开人为的干扰，可以将地震波雷达的运行时间段安排在深夜里.

图 7

Fig. 7

图 7 监测点信号的互相关系数Fig. 7 Cross-correlation coefficients of the monitoring sites′ signals

用短时傅里叶分析全部检测点未叠加的数据，在5～55 km、70 km、80～90 km和100 km处的检测点都看见了有线性调频直线，时间比地震波雷达发射信号的时间滞后(如图 8和图 9所示).在5 km检测点记录到信号里还有地震波雷达发出的二次和四次谐波信号.

图 8

Fig. 8

图 8 5 km处检测点信号图 (a)时间-振幅；(b)时间-频率.Fig. 8 Signals of the detecting site at 5 km (a)Time-Amplitude;(b)Time-Frequency.

图 9

Fig. 9

图 9 100 km处检测点信号图 (a)时间-振幅；(b)时间-频率.Fig. 9 Signals of the detecting site at 100 km (a)Time-Amplitude;(b)Time-Frequency.

在实验期间，发生了汶川地震，因此数据记录到了主震和许多余震信息，这对信号叠加极其不利，在5～55 km检测点的数据和监测点的数据，信噪比相对好一些，只采用了单次地震波雷达上升调频的数据；70 km、80～90 km和100 km检测点的数据，只进行了几次叠加.105～150 km检测点的数据，受地震信息的干扰太多，叠加后的信噪比太低，没有进行处理.重新进行数据采样，采样率由200 sps 降到50 sps.

对记录信号进行WVD变换后，再进行Hough变换，在三维极坐标内提取出角度为θ的信号，由于只关注地震波的走时和震相，对提取出的信号能量进行归一化处理，把近场监测点的波峰对应的时间作为地震波雷达发射信号的开始时间0 s(如图 10所示).

图 10

Fig. 10

图 10 监测点波的走时Fig. 10 Travel time of the monitoring site wave

再计算出每个远场检测点的地震波到达时间，扣除地震波雷达发射信号开始的参考时间后就是每个监测点地震波的走时.

中国地震局物探中心在张北地区用爆破方式布设了H-21剖面测线，该条测线与本次实验的测线相近，但并不完全重合(如图 3所示).H-21剖面测线得到了该测线地下的二维速度结构(赖晓玲等，2006)，用靠近本次实验测点的地下二维速度结构(如图 11所示)计算了折合走时曲线(Zelt and Smith， 1992；Song and Ten， 2004)(如图 12所示).H-21测线得到的直达波Pg和莫霍面反射波PmP的速度结构，用V_P/1.732和V_PmP/1.732得出直达横波V_S和莫霍面反射横波V_SmS的值.由于沽源—张北地下结构复杂，康拉德界面不连续(赖晓玲等，2006)，所以本文未计算它的直达波、反射波和折射波波速.

图 11

Fig. 11

图 11 H-21剖面图的速度结构Fig. 11 Velocity structure of H-21 section

图 12

Fig. 12

图 12 H-21剖面的折合走时Fig. 12 Reduced travel time from H-21 Section

将H-21剖面速度结构正演得到的走时曲线与本文用能量累积法提取出来的走时曲线相比较(如图 13所示)，可以看出直达纵波Pg在0～60 km范围内非常清晰，走时与用匹配滤波处理该测线的震相走时(王洪体等，2009)和H-21剖面结果完全一致；直达横波Sg在0～60 km范围内也非常清晰，走时也与用匹配滤波处理该测线的震相走时一致；莫霍面的反射纵波PmP在80 km左右开始逐渐可辨别，走时也与用匹配滤波处理该测线的震相走时和H-21剖面的走时结果一致.

图 13

Fig. 13

图 13 主要震相Fig. 13 Main seismic phases

本次地震波雷达运行控制精确，具有高度重复性，发射信号的时间和频率完全满足设计要求.用单次一个上升期的数据或简单叠加数据提取出了0～100 km内16个检测点地震波雷达信号的走时，并进行震相识别，结果表明用能量累积法提取的信号 走时与用匹配滤波法和用爆破剖面得到的结果一致，说明用能量累积法来提取地震波雷达信号走时是可靠的，用能量累积法提取出来的走时曲线更容易识别震相.

由于地下介质吸收能量，100 km以上检测点的信号相对较弱，用临时检测仪器采集记录数据，场地噪声高，再加上受汶川地震的影响，本次实验可用的数据非常有限，叠加后的信号干扰多，因此未能得到100 km以外的信号走时.本次实验测线的远端检测点和H-21测线的远端检测点不重合，而且相距越来越远，对应的地下速度结构有些差异，因此远距离检测点的震相符合程度稍差一些.结合广东新丰江的实验结果，输出力为40 t的地震波雷达信号传输有效距离在200 km左右.地震波雷达是向地下发射信号，无水平方向性，如果在方圆200 km内布置密集的高灵敏度地震仪器都能接收到信号，长期观测可以得到地壳的精细速度结构以及地震波速的动态变化情况.地震波雷达可用在城市活断层探测，水库大坝蓄水前后库区波速变化监测以及地震危险区等监测领域.

另外，本文在数据处理的过程中，未考虑信号的频散因素.在近监测点走时的主瓣和旁瓣宽带在1 s以内，提高发射信号的频带范围让旁瓣更窄，有助于提高近距离地震波走时的识别精度，这点还有待于进一步的实验.