地电阻率交流观测中信号检测方法研究*

马小溪 张 宇， 王兰炜 张兴国

1）中国北京 100085 中国地震局地壳应力研究所

2）中国北京 100085 北京市地震观测工程技术研究中心

引言

目前，我国地电前兆台网广泛使用的地电阻率观测方法是直流观测方法，其观测装置采用对称四电极装置（Zhaoetal，2011）．随着我国国民经济和社会的发展，地电阻率观测受到越来越多的电磁干扰，特别是城市轨道交通（地铁、轻轨）的干扰，其干扰频段主要集中在直流附近，且影响范围可达数十千米，严重影响了城市附近地电阻率观测台站的观测（张宇等，2014）．张世中等（2013）在北京西集、天津青光两个台站进行场地干扰测试，发现受地铁干扰的影响，两台站每天有2／3以上的小时观测值的相对方差大于0．3%，有时甚至超过3%．

采用交流供电的地电阻率方法是消除此类干扰的一个有效方法．该方法采用低频电流作为激励源，接收低频电流经大地传输后的响应信号，可以采用频率域测量技术进行选频接收（罗维斌等，2004），以获得较高的信噪比．

20世纪80年代，中国地震局曾引进美国研制的MARK型低频交流地电仪，在河北开滦马家沟地震台开展地电阻率交流法的监测试验，其结果证明了地电阻率交流观测方法的可行性（桂燮泰等，1988；马希融，1989）．

本文将首先阐述地电阻率交流观测方法的原理；然后基于地电阻率的实际测量过程，采用数值模拟方法分析不同参数条件下，数字滤波器法、频谱分析法和相关检测法等3种信号检测方法的检测效果差异；最后根据仿真结果给出适用于地电阻率交流观测方法的信号检测方法，以期为地电阻率交流观测系统的研制提供理论参考．

1 理论基础

1．1 地电阻率交流观测原理

在0．01—1 000Hz频段，交流电场分布可近似地遵循欧姆定律，电阻率交流测量方法的原理与直流测量方法相同（桂燮泰等，1988；马希融，1989；张国民等，2001）．低频交流电法测量示意图如图1所示．

图1 地电阻率交流观测方法工作原理示意图A和B为供电电极，M和N为测量电极Fig．1 Schematic diagram of principle for AC geo-resistivity method Aand Bare current-emitting electrodes，and Mand Nare measuring electrodes

地电阻率交流测量方法采用低频电流作为激励源，通过供电电极A和B向大地供入低频电流IAB，同时通过测量电极M和N测量MN之间的同频率响应信号UMN，通过计算得出地电阻率ρ的具体表达式为（张国民等，2001；刘义国等，2010）

式中K为装置系数，其计算公式为

1．2 单频信号检测法

地电阻率交流观测方法采用低频电流源发射单一频率的正弦信号，利用高采样率数据采集系统同步测量人工供电电流和人工供电电位差．要得到高精度的地电阻率观测数据，就意味着从含有噪声的测量数据中，准确检测出已知频率信号的振幅．

目前单一频率信号振幅检测方法很多，通常采用的有数字滤波器法、频谱分析法和相关检测法等3种方法．

1．2．1 数字滤波器法

数字滤波是信号处理的一种基本而重要的技术，利用数字滤波器的频率选择特性，可把滤波器的通带设置在能够覆盖有用信号的频谱范围处，使得通带内信号振幅不变，而通带外的噪声受到大幅度衰减，从而提高信噪比，得出已知频率信号的振幅值（高晋占，2007）．

数字滤波器按照设计方法分为有限冲激响应滤波器和无限冲激响应滤波器，两种滤波器理论上都能够实现高精度的滤波效果．与无限冲激响应滤波器相比，有限冲激响应滤波器有突出的优点：第一，系统总是稳定的，第二，线性相位特性；但其也有明显的缺点，即滤波器的阶数较高，在资源有限的情况下，不易满足实时需求．

常用的有限冲激响应滤波器设计方法有窗函数法、频率抽取法和最佳一致逼近法，其中窗函数设计法最为常用．对于各种窗函数来说，在滤波器阶数相同的情况下，凯塞（Kaiser）窗函数具有更好的性能（杜勇等，2012）．因此，本文选择凯塞窗来设计带通数字滤波器进行仿真．

1．2．2 频谱分析法

频谱分析法实现某一特定频率信号振幅的检测是基于任意一个函数x（t）都可以分解为无穷多个不同频率的正弦信号之和的理论，而傅立叶变换就是这样一个工具，它能通过傅里叶分析得到信号的频谱．数字信号经过傅立叶变换后可成为一个N点的复数（N为采样点数），每一个点对应着一个频率点，这个点的模值就是该频率下振幅谱的相对值．真正的振幅值为原信号相应频率下振幅的N／2倍．因此，通过傅立叶变换，可以很容易地找到指定频率信号的振幅（万永革，2012）．

离散傅立叶变换定义如下：

式中，x（n）为一有限长序列，其长度为N（门爱东等，2009）．

1．2．3 相关检测法

相关检测法是基于信号和噪声的统计特性进行信号检测，由于确定性信号不同时刻的取值之间一般都具有较强的相关性，而干扰噪声由于随机性较强，其不同时刻的取值之间相关性一般较差，利用相关性差别可以把确定性信号与干扰信号区别开来（王兰炜等，2004；张宇等，2014）．因此，从本质上说，相关函数是两个时域信号（或空间域信号）相似性的一种度量（高晋占，2007）．

对于一个已知频率的信号，可在其接收端给出与被检测信号频率相同而相位不同的正弦信号作为参考信号，将参考信号与混有噪声的输入信号进行互相关计算，从而检测出被检测信号的振幅和相位信息．相关检测法的测量原理如图2所示．含有噪声的信号x（t）与参考信号r1（t）和r2（t）相乘，分别为y1（t）和y2（t），即

图2 相关检测法的测量原理x（t）为含有噪声的信号，r1（t）和r2（t）为与x（t）频率相同而相位不同的参考信号Fig．2 Measuring principle of correlation detection method x（t）is the signals including the noise，r1（t）and r2（t）are the reference signals whose frequency are same as that of x（t）and whose phase are different from that of x（t）

式中，x（t）为被检测信号，n（t）为噪声信号，r1（t）和r2（t）为参考信号，A为被测信号振幅，ω0为被测信号频率，φ为相对相位．上式中由于信号与噪声是不相关的，因此sin（ω0t＋φ）n（t）和cos（ω0t＋φ）n（t）两项的结果为零，令这样，Y1（t）与Y2（t）经过低通滤波器后，频率为2ω0的高频项被滤波处理，可以得到z1＝AB／2cosφ和z2＝AB／2sinφ两个结果项，被检测信号的振幅值则为

2 数值模拟

2．1 输入信号的生成

采用高斯白噪声信号作为数值模拟的噪声信号，为保证模拟仿真中生成的高斯白噪声符合理想条件（均值为0，方差为1，否则会因噪声序列的不均匀而引入误差），在Matlab软件下生成长度为10万点的固定序列的高斯白噪声n（t），作为仿真的噪声序列．图3所示为10万点高斯白噪声序列，其均值为0．007 1，方差为1．005 1．

图3 10万点高斯白噪声序列（均值为0．007 1，方差为1．005 1）Fig．3 100 000-point Gaussian white noise with mean 0．007 1and variance 1．005 1

采用标准正弦信号s（t）作为被检测信号，s（t）＝Asin（2πf0t＋φ0），式中f0为信号频率，φ0为初始相位，A为被检测信号振幅，正弦信号平均功率值Ps＝A2／2．根据信噪比计算公式SNR＝10lg（Ps／Pn）保持噪声信号不变，由不同信噪比得到不同正弦信号振幅A（表1）．

输入信号x（t）由标准正弦信号和噪声两部分组成，即x（t）＝s（t）＋n（t）．图4给出了信噪比为0dB、正弦信号频率为1Hz的输入信号x（t）的波形．

2．2 数值模拟方法说明

数值模拟时，使用系统相对误差ε和相对均方差δ作为评判标准．理论上，相对误差ε越小，说明检测结果与真值越接近，准确度越高；相对均方差δ越小，说明检测结果精密度越好．因此，应该综合考虑相对误差和相对均方差结果，从数字滤波器法、相关检测法和频谱分析法中选取相对误差和相对均方差结果均较小的方法作为最优方法．

数值模拟时，选取信噪比、采样率和采样点数作为仿真参数，使用单一变量法改变参数数值，得到不同条件下的相对误差和相对均方差结果．其中，采样点数为N时，将10万点数据平均分为n＝100 000／N段，分别对每段序列进行信号检测，将得到的n个结果进行平均，其结果作为检测振幅A′＝，差值ΔA＝A′－A．

表1 信噪比与正弦信号振幅对应表Table 1 SNRs corresponding to sine signal amplitudes

均方差σ的计算公式为

相对误差ε计算公式为

相对均方差结果为均方差与检测均值的比值，计算公式为

数值模拟时，数字滤波器法仿真中设计的带通滤波器性能如表2所示．

图4 信噪比为0dB、信号频率为1Hz的输入信号x（t）的波形Fig．4 Signal with Gaussian white noise（SNR is 0dB，signal frequency is 1Hz）

表2 本文仿真中使用的带通滤波器的性能参数Table 2 Bandpass filter performance parameters used in simulation

3 仿真结果

3．1 不同信噪比的仿真结果比较

通过模拟仿真给出一定采样率和采样点数（采样时间）的条件下，不同信噪比时数字滤波器法、相关检测法和频谱分析法的检测效果．

假定采样率为10Hz，采样点数为2 000点，在信噪比为－10—50dB的情况下，表3—5分别给出了信号频率为0．1，0．5和1Hz时的仿真结果．图5给出了不同频率时的相对误差和相对均方差随信噪比变化的结果．可以看出，随着信噪比的增大，3种信号检测方法得到的相对误差和相对均方差均呈下降趋势，表明信噪比越大，信号检测效果越好．当信噪比大于20dB时，相对误差和相对均方差均趋于平稳，且不同方法的检测结果相近；当信噪比小于20dB时，频谱分析法的相对误差和相对均方差结果最小，数字滤波器法的结果最大．

从表3—5及图5中可以看出，当信号频率一定时，相关检测法与频谱分析法的差值和均方差均不随信噪比的变化而变化，说明当噪声序列固定时，相关检测法和频谱分析法的检测结果仅受到噪声序列的影响，与正弦信号的振幅大小无关；当信号频率一定时，数字滤波器法的检测差值随信噪比的增大逐渐变小，这是因为数字滤波器法采用带通滤波器对输入信号进行滤波．由于数字滤波器阶数的限制，数值模拟中滤波器的通带波动不能够达到理想状态，因此在不同频点处数字滤波器对输入信号会存在不同程度的抑制作用．根据本文仿真中选用的滤波器，其通带波动为0．01dB，说明滤波器对不同频点的抑制在0—0．01dB之间波动．假设输入信号频率为f0，振幅为A0，在频点f0处噪声的影响幅度为B，滤波器的抑制为δ，那么最终数字滤波器法检测出的信号振幅为

表3 采样率为10Hz，采样点数为2 000，信号频率为0．1Hz时的仿真结果Table 3 Simulation results of three signal detection methods on the condition of different SNR（sampling rate is 10Hz，number of samples is 2 000，signal frequency is 0．1Hz）

表4 采样率为10Hz，采样点数为2 000，信号频率为0．5Hz时的仿真结果Table 4 Simulation results of three signal detection methods on the condition of different SNR（sampling rate is 10Hz，number of samples is 2 000，signal frequency is 0．5Hz）

表5 采样率为10Hz，采样点数为2 000，信号频率为1Hz时的仿真结果Table 5 Simulation results of three signal detection methods on the condition of different SNR（sampling rate is 10Hz，number of samples is 2 000，signal frequency is 1Hz）

图5 3种方法在不同频率条件下的相对误差（左）和相对均方差（右）随信噪比的变化结果比较（采样率为10Hz，采样点数为2 000）Fig．5 Comparison of relative error（left）and relative standard deviation（right）with the variation of SNR on the condition of different frequency for three signal detection methods（sampling rate is 10Hz，number of samples is 2 000）

那么，检测得到的误差结果为

因此，当信噪比变化，即被检测信号振幅变化时，由于所使用的带通滤波器没有变，通带范围内噪声的影响幅度不变，滤波器对已知频点的抑制作用不变，由式（9）可知，数字滤波器法的检测误差将随信噪比的增大而减小．

3．2 不同采样率的仿真结果比较

通过模拟仿真，给出了在一定采样时间和信噪比的条件下，当采样率不同时，数字滤波器法、相关检测法和频谱分析法的检测效果．

假定采样时间为200s，信噪比为－10dB，采样率为10—30Hz，表6和表7分别给出了信号频率为0．1Hz和1Hz时的相对误差和相对均方差的仿真结果．可以看出，在采样率为10—30Hz、信号频率分别为0．1Hz和1Hz时，频谱分析法的相对误差结果和相对均方差结果均最小，信号检测效果要优于相关检测法和数字滤波器法．

3．3 不同采样点数的仿真结果比较

通过模拟仿真，给出了在一定的采样率和信噪比条件下，当采样点数（采样时间）和信号频率不同时，数字滤波器法、相关检测法和频谱分析法的检测效果．

假定采样率为10Hz，信噪比为－10dB，采样点数为2 000—10 000点，表8和表9分别给出了信号频率为0．1Hz和1Hz时的相对误差和相对均方差的仿真结果．可以看出：随着采样点数的增加，相关检测法和频谱分析法的差值变化不大，频谱分析法的相对误差结果要优于其它两种方法；随着采样点数的增加，信号检测的均方差结果逐渐减小，其原因是采样点数的增加，能够提高频谱的分辨能力．

表7 采样时间为200s，信噪比为－10dB，信号频率为1Hz时的仿真结果Table 7 Simulation results of three signal detection methods on the condition of different sampling rate（SNR is－10dB，sampling time is 200s，signal frequency is 1Hz）

表8 采样率为10Hz，信噪比为－10dB，信号频率为0．1Hz时的仿真结果Table 8 Simulation results of three signal detection methods on the condition of different sampling number（SNR is－10dB，sampling rate is 10Hz，signal frequency is 0．1Hz）

表9 采样率为10Hz，信噪比为－10dB、信号频率为1Hz时的仿真结果Table 9 Simulation results of three signal detection methods on the condition of different sampling number（SNR is－10dB，sampling rate is 10Hz，signal frequency is 1Hz）

从上述仿真结果中可以得出以下结论：

1）相同条件下，信噪比越高，信号检测的相对误差和相对均方差越小，检测效果越好．信噪比为－10—50dB时，频率分析法和相关检测法的信号检测结果要优于数字滤波器法；当信噪比高于20dB时，3种信号检测方法的检测效果差别不大，且随着信噪比的提高，相对误差和相对均方差的变化趋于平稳．

2）相同条件下，采样率为10—30Hz时，频谱分析法的相对误差和相对均方差最小，信号检测效果优于相关检测法和数字滤波器法．

3）相同条件下，随着采样点数的增加，信号检测的均方差逐渐减小，相关检测法和频谱分析法的差值变化不大．从相对误差结果上看，频谱分析法和相关检测法的检测效果要优于数字滤波器法．

4）信号频率一定时，对于同一噪声序列，相关检测法和频谱分析法的检测结果仅与噪声序列的频谱相关，与正弦信号的振幅无关．

4 讨论与结论

本文利用数值方法模拟频谱分析法、相关检测法和数字滤波器法等3种信号检测方法，得出在不同参数条件下对单一频率输入信号的检测精度．通过对比分析检测结果，并结合3种方法在实际系统中实现的难易程度，给出适用于地电阻率交流观测系统的最优信号检测方法．

从仿真结果上看，频谱分析法和相关检测法的相对误差和相对均方差结果均优于数字滤波器法．其主要原因在于，数字滤波器法中所使用的带通滤波器有一定的通带范围，所受到的噪声影响较大；另外，受滤波器阶数和仿真参数的限制，带通滤波器的纹波因数不能达到理想状态，对被测信号会产生一定程度的抑制．两种影响的同时作用使得数字滤波器法的检测效果不够理想．

相关检测法的检测效果与频谱分析法差别不大．但是，相关检测法在使用时需要满足两个条件：① 该检测方法中的参考信号需要与被测信号的频率完全一致，这就要求在实际应用时低频电源发射的供电信号需要非常稳定；② 该检测方法中需要进行低通滤波，低通滤波器的通带范围应该在直流附近，通带越窄越好．若采用有限冲激响应滤波器设计低通滤波器的话，会需要很高阶数的滤波器，从而导致相关检测法需要占用更多的系统资源．并且由于滤波器性能的限制，相关检测法的检测效果也会受到影响．因此，在实际应用中，该检测方法的实现难度较大，占用的系统资源也较多．

与相关检测法相比，频谱分析法的处理方法更简单，易于实现．但是，该检测方法的前提是噪声频率与被提取信号的频率不在同一个频段．因此，在实际观测时，需要提前对场地进行测试，找到干扰比较小的频段，选择合适的信号频率．

目前直流地电阻率观测的相对误差要求不大于0．3%（中国地震局，2009），因此采用地电阻率交流测量方法时，在采样率不变的前提下，通过提高采样时间长度，以增加频谱分析时的数据长度，提高频率分辨力；或者采用滑动多次平均的方法来提高频谱分析法对数据处理的准确度，以提高检测精度．

综上所述，频谱分析法是适用于地电阻率交流观测的一种较好的信号检测技术，希望能够为地电阻率交流观测系统的研制提供借鉴．

杜勇，路建功，李元洲．2012．数字滤波器的MATLAB与FPGA实现［M］．北京：电子工业出版社：120-130．

Du Y，Lu J G，Li Y Z．2012．ImplementationofDigitalFilterBasedonMATLABandFPGA［M］．Beijing：Publishing House of Electronics Industry：120-130（in Chinese）．

高晋占．2007．微弱信号检测［M］．北京：清华大学出版社：3，239-279．

Gao J Z．207．DetectionofWeakSignals［M］．Beijing：Tsinghua University Press：3，239-279（in Chinese）．

桂燮泰，戴经安，关华平．1988．低频交流电法的试验和研究［J］．西北地震学报，10（2）：22-28．

Gui X T，Dai J A，Guan H P．1988．Experiment and discussion on the low frequency AC method［J］．NorthwesternSeismologicalJournal，10（2）：22-28（in Chinese）．

刘义国，董浩斌，谭超．2010．基于LabVIEW的低频交流电法仪的设计［J］．工矿自动化，36（5）：88-91．

Liu Y G，Dong H B，Tan C．2010．Design of low-frequency AC method instrument based on LabVIEW［J］．Industry andMineAutomation，36（5）：88-91（in Chinese）．

罗维斌，白宜诚，杨学顺．2004．用交流电阻率法探测煤矿导（含）水构造［J］．物探与化探，28（2）：139-141．

Luo W B，Bai Y C，Yang X S．2004．The application of the AC resistivity method to the exploration of water-bearing and water-conducting structures in coal mines［J］．Geophysical＆GeochemicalExploration，28（2）：139-141（in Chinese）．

马希融．1989．交流地电阻率方法在马家沟地震台试验结果［J］．地震，（3）：46-51．

Ma X R．1989．Experimental results of alternative resistivity method in Majiagou station［J］．Earthquake，（3）：46-51（in Chinese）．

门爱东，苏菲，王雷，王海婴，李江军．2009．数字信号处理［M］．第二版．北京：科学出版社：75-136．

Men A D，Su F，Wang L，Wang H Y，Li J J．2009．DigitalSignalProcessing［M］．2nd ed．Beijing：Science Press：75-136（in Chinese）．

万永革．2012．数字信号处理的MATLAB实现［M］．第二版．北京：科学出版社：86-132．

Wan Y G．2012．ImplementationofDigitalSignalProcessingBasedonMATLAB［M］．2nd ed．Beijing：Science Press：86-132（in Chinese）．

王兰炜，赵家骝，王子影，王燕琼．2004．相关检测在甚低频电磁信号检测中的应用［J］．西北地震学报，26（4）：339-342．

Wang L W，Zhao J L，Wang Z Y，Wang Y Q．2004．Application of correlation detection to ELF signal detection［J］．NorthwesternSeismologicalJournal，26（4）：339-342（in Chinese）．

张宇，王兰炜，张兴国，朱旭，刘大鹏，颜蕊．2014．相关检测技术在低频交流地电阻率观测中的应用［J］．地球物理进展，29（4）：1973-1979．

Zhang Y，Wang L W，Zhang X G，Zhu X，Liu D P，Yan R．2014．Application of correlation detection technology in lowfrequency AC geo-resistivity observation［J］．ProgressinGeophysics，29（4）：1973-1979（in Chinese）．

张国民，傅征祥，桂燮泰．2001．地震预报引论［M］．北京：科学出版社：241-246．

Zhang G M，Fu Z X，Gui X T．2001．IntroductiontoEarthquakePrediction［M］．Beijing：Science Press：241-246（in Chinese）．

张世中，石航，王兰炜，胡哲，刘大鹏，魏连生，鞠永．2013．地电台站受城市轨道交通干扰的测试分析与抗干扰措施研究［J］．地震学报，35（1）：117-124．

Zhang S Z，Shi H，Wang L W，Hu Z，Liu D P，Wei L S，Ju Y．2013．Test analysis on disturbances caused by urban rail transit at geoelectric stations and measures to reduce its influence［J］．ActaSeismologicaSinica，35（1）：117-124．

中国地震局．2009．地震标准汇编［M］．北京：地震出版社：1214-1220．

China Earthquake Administration．2009．SeismicStandardsCompilation［M］．Beijing：Seismological Press：1214-1220（in Chinese）．

Zhao J L，Wang L W，Qian J D．2011．Research on geo-electrical resistivity observation system specially used for earthquake monitoring in China［J］．EarthquakeScience，24（6）：497-511．