地球物理学报  2015, Vol. 58 Issue (4): 1259-1268   PDF    
引用本文
刘明辉, 周银兴, 彭朝勇, 陈阳, 李江, 王洪体. 2015. 用能量累积法检测地震波雷达信号. 地球物理学报, 58(4): 1259-1268, doi: 10.6038/cjg20150414   
LIU Ming-Hui, ZHOU Yin-Xing, PENG Chao-Yong, CHEN Yang, LI Jiang, WANG Hong-Ti. 2015. Using the method of accumulating energy for detecting seismic wave radar signal. Chinese J. Geophysics, 58(4): 1259-1268, doi: 10.6038/cjg20150414   
用能量累积法检测地震波雷达信号
刘明辉1,2, 周银兴1,2, 彭朝勇3, 陈阳1,2, 李江1,2, 王洪体1,2    
1. 中国地震局地震预测重点实验室, 北京 100036;
2. 中国地震局地震预测研究所, 北京 100036;
3. 中国地震局地球物理研究所, 北京 100081
投稿时间:2013-12-13    最后修改时间: 2015-03-30 .
基金项目：国家科技支撑计划专题(2012BAF14B12-04)、国家自然科学基金项目(41404048)、中国地震局地震预测研究所基本科研业务专项(2010IES0201)和中国地震局地球物理研究所基本科研业务专项(DQJB14B05)联合资助.
作者简介：刘明辉,男,1976年生,工程师,现从事地震观测技术的研究、推广和应用,E-mail:geodevice@sina.com
通讯作者：彭朝勇,男,1982年生,助理研究员,现从事地震应急观测技术研究. E-mail:pengchaoyong@cea-igp.ac.cn
摘要：将具有高度重复性的地震波雷达长时间向地壳内发射线性调频信号,经过地下介质的传播后到达地面,用地震仪器在监测点和检测点记录下来,通过分析数据来了解地壳速度结构及波速变化.线性调频信号是一种非平稳信号,它的频率随时间线性变化,有很好的能量聚集性,非常适合做时间-频率分析.本文用短时傅里叶变换对监测点的信号进行时间-频率分析,以检验地震波雷达发射信号的时间和频率是否和控制系统一致.通过Wigner-Ville分布将地震波雷达发射的信号能量聚集在线性调频直线上,再用Hough变换累积聚集的能量形成波峰,按照线性调频直线的倾角提取波峰所在行,计算到时后构成地震波走时曲线图.用靠近本次实验地点的H-21剖面得到的地壳速度结构正演该测线的Pg、Sg、PmP和SmS的折合走时曲线,并与用能量累积法提取出的地震波走时曲线进行对比,分析结果表明:地震波雷达发射线性调频信号的时间和频率都符合控制要求,重复性高达99.9%以上,可以清晰地分辨出Pg、Sg震相,并且PmP和SmS震相可辨.
关键词： 地震波雷达     线性调频     能量聚集     能量累积     波速结构     波速变化     绿色环保    
Using the method of accumulating energy for detecting seismic wave radar signal
LIU Ming-Hui1,2, ZHOU Yin-Xing1,2, PENG Chao-Yong3, CHEN Yang1,2, LI Jiang1,2, WANG Hong-Ti1,2    
1. Key Laboratory of Earthquake Prediction, China Earthquake Administration, Beijing 100036, China;
2. Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China;
3. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
Abstract: How to identify seismic wave velocity anomalies and to study their characteristics has practical significances for earthquake study and prediction. Almost all disastrous earthquakes in China's mainland have occurred in the Earth's crust within depth of 5~25 km. Therefore, dynamic monitoring physical parameters of the Earth's crust is an effective approach to predicting earthquakes. Seismic Wave Radar (SWR) is a type of mechanical and electrical devices. It continuously excites Linear Frequency Modulation (LFM) signals into the Earth's crust, and these signals are recorded by high sensitivity seismographs at the deployment sites. How to retrieve impulsive seismic waveforms from the recorded data which include signals and different kinds of noises is very important for calculating accurate travel times and temporal velocity variation. In this work, we propose a new method for processing the SWR data, which is based on accumulating energy in time-frequency domain. LFM is a non-stationary signal and its frequencies vary linearly as a function of time and have a good energy concentration, so it is very suitable for time-frequency analysis. We use Short-Time Fourier Transform (STFT) to analyze the data recorded by the monitoring stations, and acquire the time-frequency distribution, in order to validate whether the excitation time and frequencies are consistent with those set by the control system of SWR. Then, two methods are introduced for the waveform retrieval from the SWR data. One is the Wigner-Ville Distribution (WVD) algorithm with the best time-frequency concentration capability for the LFM signals, and the other one is the Wigner-Hough Transform (WHT) which is helpful to suppress the cross-time interference in the signal detection and parameter estimation for the multi-component LFM signals. By aggregating energy of the excited signals of the SWR on the LFM line with WVD method and accumulating the concentration energy to a crest by WHT, we can extract the peak row according to the inclination angle of the LFM line and compose the seismic wave travel-time curve after calculating the travel time.Firstly, in order to quantify the repeatability of the SWR, we calculated cross-correlation coefficients of all waveform pairs retrieved from different time windows. Within the total 126 hours' SWR data, cross-correlations coefficients of 123 hours' are bigger than 0.999. The high correlations between the waveforms indicate the excellent repeatability of the SWR. In addition, the results obtained from STFT methods show that the excitation time and frequencies of the LFM signals met the control requirements and the repeatability is more than 99.9%. Finally, the crustal velocity structure obtained from the H-21 section near the location of this experiment was used to calculate the reduced travel-time curves of Pg, Sg, PmP and SmS of the detecting line, which were compared with those retrieved by our method. The results demonstrate that the proposed method is an effective way for retrieving waveform, identifying seismic phase and measuring travel time. By using this method, we can clearly identify Pg and Sg, and high frequency seismic phases such as PmP and SmS with strong amplitude in some epicenter distances can be discernable. In all, the SWR with high repeatability can be easily applied to monitoring temporal changes and imaging the spatial variations of the subsurface structure. The proposed method provides a feasible way to process the SWR data and advance its application in monitoring the crustal processes.
Key words: SWR     LFM     Energy concentration     Accumulating energy     Wave velocity structure     Wave velocity variation     Green    
1 引言

当地震波通过应力变化区时其传播速度也要发生变化，由于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 m2，相对天然地震而言，它是一个很准确的“点”源.地震波雷达是机电设备，控制精确、功耗低、噪音小、重复性高、无破坏性、可移动使用，是一种绿色环保的震源.俄罗斯在西伯利亚，距地震波雷达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).

2 数据 2.1 地震波雷达信号

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

图 1 地震波雷达图
(a)地震波雷达外观结构；(b)地震波雷达运行示意图.Fig. 1 Seismic wave radar
(a)Appearance of the SWR structure;(b)SWR running diagram.
式中，F(t)为向地面输出的垂直力，M为偏心轮的质量，r为偏心轮的转动半径，ω(t)为偏心轮转动的角频率，φ(t)为偏心轮转动相对于固定位置的相位.M，r 是地震波雷达的机械参数，在仪器制造时就已经 测得，因此它输出力的大小和相位都取决于 ω(t).

假设线性调频的起始频率为fstart，结束频率为fend，调频周期为T，单位为s，则调频斜率为fk：

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

2.2 发射信号的调频模式

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

图 2 地震波雷达的调频模式Fig. 2 Frequency modulation model of SWR
2.3 监测数据与检测数据

项目期间，在河北沽源县安装了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 监测点和检测点分布图(王洪体等，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.

3 信号分析方法

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

3.1 短时傅里叶变换

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

式中，h(t)为窗函数.一般情况下，h(t)为紧支集的实函数且其傅里叶变换的能量集中在低频处，可以看作低通滤波器的脉冲响应.从式(3)可以看出，在短时傅里叶变换中，首先利用平移参数t，通过窗函数h(t)对信号s(t)进行分段，即

获得原信号在t时刻附近τ时段的信号，然后对该段信号st(τ)进行傅里叶变换.随着平移参数t的改变，h(t)所确定的“时间窗”在时间轴上移动，使得信号s(t)“逐步”进入被分析的状态.由于分段信号st(τ)突出了原信号在t时刻附近一个时段的特性，所以分段信号st(τ)的傅里叶变换反映了原信号在该时段的频率分布.

3.2 聚集信号能量

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

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

图 4 监测点信号的WVDFig. 4 WVD of the monitoring site signal
3.3 累积信号能量

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

Hough变换把直角坐标表示的图像中的每一个点(x，y)变换为极坐标中的一条相关的正弦曲线，曲线上每一点的幅度等于图像中象素点(x，y)的亮度.图像中的所有点通过Hough变换转换为极坐标中的相交的一束正弦曲线，也就是说，Hough变换在图像中沿着直线积分，每个积分值都会影响与这条直线的参数相关联的点(ρ，θ).如果图像中一些高亮度值的像素点聚集在一条直线附近，那么将在极坐标平面内观察到一个波峰，其坐标与直线的参数直接相关(胡昌华等，2002).

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

图 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变换后，在三维极坐标内角度为θ的信号提取出来，把近场监测点波峰对应的时间作为地震波雷达发射信号开始的参考时间，计算出每个远场检测点的线性调频信号到达时间，这就是通过能量累计法来实现对地震波雷达信号的检测方法.

4 结果分析 4.1 监控点信号的时间-频率分析

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

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

为了监测地下介质波速的变化情况，要用地震波雷达长期向地下发射信号，地震波雷达运行的重复性显得尤为重要.信号相关为波形相似的度量，通常用相关系数来表示，相关系数越大，则表明波形就越相似，如果两个信号的相关系数为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 Cross-correlation coefficients of the monitoring sites′ signals
4.3 检测点信号时间-频率分析

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

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

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

在实验期间，发生了汶川地震，因此数据记录到了主震和许多余震信息，这对信号叠加极其不利，在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 Travel time of the monitoring site wave

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

4.5 震相识别

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

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

图 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 Main seismic phases
5 结论

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

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

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

致谢 感谢中国地震局地震预测研究所庄灿涛研究员对本文研究成果的指导.感谢北京港震机电技术有限公司承担了地震波雷达的运行维护和监测任务以及为本文的研究提供的数据支持.感谢评审人提出的宝贵意见，使得本文的研究内容得到了进一步的提高.

参考文献
[1] Alekseev A S, Chichinin I S, Korneev V A. 2005. Powerful low-frequency vibrators for active seismology. Bulletin of the Seismological Society of America, 95(1):1-17.
[2] Chen Y, Zhang W, Chen H L, et al. 2006. Seismic radar. Progress in Geophys. (in Chinese), 21(1):1-5.
[3] Cheng Z L, Zheng S H, Liu J. 2007. Application Research of Digital Seismological Observation Data (in Chinese). Beijing:Seismological Press.
[4] Department of science & earthquake monitoring of China Earthquake Administration. 1995. Earthquake Observation Technology (in Chinese). Beijing:Seismological Press.
[5] Feng R, Pang Q Y, Fu Z X, et al. 1976. Variations of Vp/Vs before and after the Haicheng Earthquake of 1975. Chinese J. Geophys. (in Chinese), 19(4):295-305.
[6] Feng R. 1977. On the variations of the velocity ratio before and after the Xinfengjiang reservoir impounding earthquake of M=6.1. Chinese J. Geophys. (in Chinese), 20(3):211-221.
[7] Hu C H, Zhou T, Xia Q B, et al. 2002. System Analysis & Design Based on MATLAB-Time-Frequency (in Chinese). Xi'an:Xi'an Electronic and Science University Press.
[8] Ikuta R, Yamaoka K. 2004. Temporal variation in the shear wave anisotropy detected using the Accurately Controlled Routinely Operated Signal System (ACROSS). Journal of Geophysical Research, 109(B9):B09305, doi:10.1029/2003JB002901.
[9] Lai X L, Zhang X K, Sun Y. 2006. The complexity feature of crust-mantle boundary in Zhangbei seismic region and its tectonic implication. Acta Seismologica Sinica (in Chinese), 28(3):230-237.
[10] Li M X, Liu J. 2006. Study on velocity ratio (Vp/Vs) anomaly of earthquake sequences in Yunnan region. Earthquake (in Chinese), 26(1):26-34.
[11] Li Z W, You Q Y, Ni S D, et al. 2013. Waveform retrieval and phase identification for seismic data from the CASS experiment. Pure and Applied Geophysics, 170(5):815-830.
[12] Liu X K, Cui R S, Wang H T, et al. 2013. Detecting of controlled accurate seismic source signal using Wigner-Hough transformation. Earthquake (in Chinese), 33(3):33-42.
[13] Peng C Y, Yang J S, Xue B, et al. 2013. Research on correlation between early-warning parameters and magnitude for the Wenchuan Earthquake and its aftershocks. Chinese J. Geophys. (in Chinese), 56(10):3404-3415, doi:10.6038/cjg20131016.
[14] Peng C Y, Zhu X Y, Yang J S, et al. 2013a. Development of an integrated onsite earthquake early warning system and test deployment in Zhaotong, China. Computers & Geosciences, 56:170-177, doi:10.1016/j.cageo.2013.03.018.
[15] Peng C Y, Yang J S, Xue B, et al. 2014a. Exploring the feasibility of earthquake early warning using records of the 2008 Wenchuan earthquake and its aftershocks. Soil Dynamics and Earthquake Engineering, 57:86-93, doi:10.1016/j.soildyn.2013.11.005.
[16] Peng C Y, Yang J S, Zheng Y, et al. 2014b. A τc magnitude estimation of the 20 April 2013 Lushan earthquake, Sichuan, China. Science China Earth Sciences, 57(12):3118-3124, doi:10.1007/s11430-014-4971-8.
[17] Peng C Y, Yang J S, Xue B, et al. 2015. A low-latency data acquisition system for earthquake early warning. Earthquake (in Chinese), 35(1):140-148.
[18] Song J L, Ten Brink U. 2004. RayGUI 2.0:A Graphical User Interface for Interactive Forward and Inversion Ray-Tracing:USGS Open-File Report 2004-1426. New York:BiblioGov.
[19] Tang X H, Li Q L. 2008. Time-frequency & wavelet transfer (in Chinese). Beijing:Science Press.
[20] Wang H T, Zhuang C T, Xue B, et al. 2009. Precisely and actively seismic monitoring. Chinese J. Geophys. (in Chinese), 52(7):1808-1815, doi:10.3969/j.issn.0001-5733.2009.07.015.
[21] Yang W, Ge H K, Wang B S, et al. 2010. Velocity changes observed by the precisely controlled active source for the Mianzhu Ms5.6 earthquake. Chinese J. Geophys. (in Chinese), 53(5):1149-1157, doi:10.3969/j.issn.0001-5733.2010.05.016.
[22] Yang W, Wang B S, Ge H K, et al. 2013. Characteristics and signal detection method of accurately controlled routinely operated signal system. Journal of China University of Petroleum (Edition of Natural Science) (in Chinese), 37(1):50-55, 69, doi:10.3969/j.issn.1673-5005.2013.01.008.
[23] Zelt C A, Smith R B. 1992. Seismic traveltime inversion for 2-D crustal velocity structure. Geophys. J. Int., 108(1):16-34.
[24] Zhang X D, Bao Z. 1998. Non staionary signal analysis & processing (in Chinese). Beijing:National Defense Industry Press.
[25] Zhao C L, Lu C, Hao T Y, et al. 2013. A study of the high-precision modular lightweight small vibrator. Chinese J. Geophys. (in Chinese), 56(11):3690-3698, doi:10.6038/cjg20131110.
[26] Zhuang C T. 2002. A technical measure to study and predict disastrous earthquakes occurring directly beneath big cities. Recent Developments in World Seismology (in Chinese), (8):35-37.
[27] 陈颙, 张蔚, 陈翰森等. 2006. 地震雷达. 地球物理学进展, 21(1): 1-5.
[28] 陈章立, 郑斯华, 刘杰. 2007. 数字地震观测资料应用研究. 北京: 地震出版社.
[29] 冯锐, 庞庆衍, 傅征祥等. 1976. 海城地震前后地震波速比的变化. 地球物理学报, 19(4): 295-305.
[30] 冯锐. 1977. 新丰江6.1级水库地震前后的波速比变化. 地球物理学报, 20(3): 211-221.
[31] 国家地震局科技监测司. 1995. 地震地形变观测技术. 北京: 地震出版社.
[32] 胡昌华, 周涛, 夏启兵等. 2002. 基于MATLAB的系统分析与设计——时频分析. 西安: 西安电子科技大学出版社.
[33] 赖晓玲, 张先康, 孙译. 2006. 张北地震区壳幔边界复杂性特征及其构造意义. 地震学报, 28(3): 230-237.
[34] 黎明晓, 刘杰. 2006. 云南地区地震序列的波速比(Vp/Vs)异常研究. 地震, 26(1): 26-34.
[35] 刘希康, 崔仁胜, 王洪体等. 2013. 用Wigner-Hough变换检测精密可控震源信号. 地震, 33(3): 33-42.
[36] 彭朝勇, 杨建思, 薛兵等. 2013. 基于汶川主震及余震的预警参数与震级相关性研究. 地球物理学报, 56(10): 3404-3415, doi: 10.6038/cjg20131016.
[37] 彭朝勇, 杨建思, 薛兵等. 2015. 用于地震预警的通用数据采集系统构建. 地震, 35(1): 140-148.
[38] 唐向宏, 李奇良. 2008. 时频分析与小波变换. 北京: 科学出版社.
[39] 王洪体, 庄灿涛, 薛兵等. 2009. 精密主动震源监测. 地球物理学报, 52(7): 1808-1815, doi: 10.3969/j.issn.0001-5733.2009.07.015.
[40] 杨微, 葛洪魁, 王宝善等. 2010. 由精密控制人工震源观测到的绵竹5.6级地震前后波速变化. 地球物理学报, 53(5): 1149-1157, doi: 10.3969/j.issn.0001-5733.2010.05.016.
[41] 杨微, 王宝善, 葛洪魁等. 2013. 精密控制机械震源特征及信号检测方法. 中国石油大学学报(自然科学版), 37(1): 50-55, 69, doi: 10.3969/j.issn.1673-5005.2013.01.008.
[42] 张贤达, 保铮. 1998. 非平稳信号分析与处理. 北京: 国防工业出版社.
[43] 赵春蕾, 卢川, 郝天珧等. 2013. 高精度组合式轻便小型可控震源的研究. 地球物理学报, 56(11): 3690-3698, doi: 10.6038/cjg20131110.
[44] 庄灿涛. 2002. 探索预报大城市直下型灾害性地震的一种技术措施. 国际地震动态, (8): 35-37.