引言

地电阻率观测是开展地震预测预报的有效方法(赵玉林, 钱复业, 1978, 1979；钱复业，赵玉林，1980; 钱家栋，1993；钱家栋等，2013). 1966年以来，我国的地震地电阻率观测经历了从地表测量到井下观测的发展历程，在观测技术、理论方法、数据应用等方面均取得了较大进展(赵和云，1994; 赵家骝等，1996；汤吉等, 1998；Ma, 2002；Lu et al, 2004, 2016；黄清华，2005；Huang, 2008；张学民等，2009；杜学彬, 2010；高曙德等, 2010, 2013, 2016；王兰炜等, 2011；Zhao et al，2011；赵国泽等，2012；杜学彬等, 2015).目前, 地震台站受到日益严重的电磁干扰，对地电阻率观测造成了巨大的影响，干扰引起的视电阻率变化淹没了震前孕育的地电阻率变化(高曙德，2016).地震台站亟需引入新的观测技术和方法来提高地电阻率观测精度.

相关辨识法是一种可以有效去除随机噪声和干扰的系统辨识方法. Duncan等(1980)最早尝试了基于伪随机二进制序列(pseudo random binary sequence，简写为PRBS，简称m序列)为信号源的宽频带地电阻率测量方法，证明采用编码源进行大深度地电阻率观测的技术路线是可行的. Ziolkowski等(2007)基于PRBS编码源技术开发出了多道瞬变电磁法，对大埋深高阻薄层有较高的探测能力.国内近几年也掀起对基于PRBS编码源电磁勘探法理论和应用的研究热情. 赵璧如等(2006)介绍了基于CDMA(code division multiple access)技术的地电阻率高精度测量方法；汤井田和罗维斌(2008)以及罗维斌等(2012, 2016)提出了基于逆重复M序列伪随机频率域电磁勘探方法原理，并讨论了对各种干扰的压制方法以及提高测量精度的措施.逆重复M序列由2个PRBS序列(m序列)隔位取反产生，周期是m序列的2倍，自相关峰正反向交替，无直流分量，所以逆重复M序号是系统辨识一种比较理想的信号源.基于PRBS伪随机电法及观测系统研究也有新的进展(淳少恒等，2014；罗先中等，2014；罗延钟等, 2015, 2016)；齐彦福等(2015)利用方波响应移位叠加和电流导数与阶跃响应褶积两种方法实现了理论m序列和实际发射波形的全时正演模拟，武欣等(2015)研究了m序列伪随机编码源电磁响应的精细辨识问题；王若等(2016)利用伪随机编码源激发下的时域电磁信号合成与理论值对比，可用于检验后续数据处理方法的可靠性和有效性；张世中等(2001)采用有限长的伪随机码来控制供电周期(变周期供电测量法)以减少测量的随机误差, 并将其应用于ZD8B地电仪提高观测系统的测量精度；张宇等(2014, 2016)利用相关检测技术研制了交流地电阻率观测仪器，并在江宁台进行观测试验，其结果显示抗干扰能力得到了改善，但频段相对较窄.

近年来编码源电阻率测量方法在地学勘探领域得到了快速发展且应用广泛，主要采用m序列为激发信号源.该方法在地震系统目前尚处于试验阶段, 主要原因是没有配套的观测仪器，其次是现阶段地震系统原有的观测线路和观测场地无法改造，该方法能否观测到可靠连续的数据，还需要作进一步试验研究.本文基于2ⁿ系列伪随机电磁法理论(何继善，2010)，采用循环互相关法获取系统频率响应的方法(相关辨识方法是用伪随机序列(m序列或逆重复M序列)作为激励信号输入系统，由观测电场响应和输入电流信号辨识得到的地电系统频率响应，反映待探测大地系统的阻抗信息)，所获取的阻抗谱可直接用于反演解释，也可反算出系统的传递函数，从而获得地电系统的阶跃响应以及冲激响应，可得到复电阻率、阻抗及Cole-Cole模型参数(罗维斌, 李梅，2015).

本文基于相关辨识方法开展仪器的研制和实用化创新，在环境干扰较大的甘肃省兰州观象台和陇南汉王地震台利用现有的观测场地和线路进行试验观测，获取相关数据，以检验该方法在强干扰环境下的抗干扰能力，以期拓宽直流电法测量范围，提高地电阻率的观测精度.

1.   编码源理论和循环互相关法

1.1   基本原理

采用PRBS或逆重复M序列作为激励信号波形，利用互相关法辨识未知系统的频率响应，系统输入、输出的互相关为(刘立超等，2014；罗维斌等，2016)：

式中，h_s(τ)为观测系统的冲激响应，h_e(τ)为待辨识未知系统的冲激响应，R_u(τ)为输入信号的自相关，R_yu(τ)为系统输入与输出信号的互相关函数.在复数域，H_s(s)和H_e(s)分别为观测系统自身h_s(τ)和待辨识系统h_e(τ)的拉普拉斯变换, 式(1) 在频率域可变换为

式中，W_yu(jω)为系统输入与输出信号的互相关函数R_yu(τ)的拉普拉斯变换，W_u(jω)为系统输入信号自相关函数R_u(τ)的拉普拉斯变换，H_e(jω)为系统的频率响应.

基于循环互辨识法的工作流程如图 1所示.在计算数据中引入了一个参考信号S(按照输入激励信号的时间序列周期和采样率触发软件程序生成的无干扰信号)，与系统输入信号U(t)和系统输出信号y_out(t)作循环互相关，由循环互相关时间序列的快速傅里叶变换得到互功率谱，由式(3) 和(4) 恢复出系统的频率响应H(ω)和相位φ(ω)，这样处理的优点是即使输入和输出两端都存在工频噪声, 也能得到有效抑制.

图  1  基于逆重复M序列系统辨识技术流程图

Figure  1.  Workflow of system identification based on inverse repeated M-sequence

下载: 全尺寸图片 幻灯片

式中，H(ω)为频率响应，P_Cys为参考信号S与输出信号y_out的互功率谱，P_Cus为参考信号S与输入信号U(t)的互功率谱，φ(ω)为相位，mod为模运算符.

对于地电阻率谱而言，当观测值是接收电极间的电压值U_MN时间序列时(MN为接收电场的测量间距)，以上计算过程可改写为

式中，P_CUMN^s(ω)为参考信号S与接收电压信号U_MN的互功率谱，P_CIAB^s(ω)为参考信号S与供电电流信号I_AB(t)的互功率谱，则为由观测电压U_MN和激发电流I循环互相关辨识得到的大地阻抗谱.由此可得地电阻率谱为

式中K为观测装置系数

式中，A和B为供电电极，M和N为测量电极.由于采用互相关技术，本测量方法对观测场地干扰有较强的压制能力.

1.2   循环互相关法

对于两个时间序列的互相关计算，采用非循环互相关法计算时，在序列移位过程中需要在序列头部和尾部补零来实现两个序列长度相等，计算结果会在相关脉冲峰两侧产生误差干扰；循环互相关法计算时，将等长的两个时间序列分别首尾环状相连，这样在计算中序列循环移位，无需补零(Roinila et al, 2009)，且互相关峰值相等，这个特性在多周期信号功率谱叠加中优势显著.

2.   试验结果与分析

2.1   观测方案

利用甘肃省地震局地电阻率台站的观测装置开展试验, 仪器选用由罗维斌等(2016)试制的测量样机，数据采样率为1000 Hz，对发送电流和接收电压多周期同步采集为时间序列，研究同步采样和非同步采样对于辨识结果的影响，通过计算观测结果来调整和优化仪器参数，从而确定最佳的观测参数.

2.2   兰州观象台试验

2016年7月30日选用兰州观象台进行试验.该台站位于黄河Ⅲ级阶地(次生黄土)，第四纪覆盖层厚10—30 m，岩性为亚砂土、粉砂、沙砾石，下伏砂岩.地下水位埋深5—30 m，含水层厚度大于5 m，含水层岩性为粉砂和沙砾石; 台站的电测曲线为Q型(李乐进，1985). 1970年5月开始观测, 布设有3道测线，其方位分别为N5°W, N40°E，N50°W，近年受站外蔬菜大棚建设和工业及城镇化建设的干扰影响增大, 测量数据波动较大(图 2)，其中2015年初至2016年初，N5°W连续10个月出现脉冲状干扰，最大干扰值可达150 Ω·m，测值曲线无法清楚地看到工程焊接等干扰，所以在绘图时对图 2作了限幅来显示工程建设和焊接的干扰. 2016年初N40°E测道供电极B附近电焊作业，资料出现振荡变化；2016年9月开始建设和工程焊接使N5°W测值振荡变化持续了8个月，给资料分析带来了不确定的因素, 资料的可用性大大降低，该台目前已规划选址搬迁.本次试验利用N40°E现有观测线路(AB=900 m，MN=300 m), 用8号镀锌铁丝悬挂电缆线作为测量及供电导线，电极为铅板，埋深2.5 m.

图  2  2015年1月—2017年6月兰州观象台直流电法观测的地电阻率曲线

Figure  2.  The georesistivity curves recorded at the Lanzhou observatory by direct current method from January 2015 to June 2017

下载: 全尺寸图片 幻灯片

2.2.1   输入信号和输出电压

为了验证信号序列长度压制干扰的能力，设计编码源信号产生参数：寄存器位数m=7，序列周期254位；晶振频率5 MHz，分频后获得5.24 Hz时钟频率激励移位寄存器，产生的逆重复M序列周期为63.5 s.直流电压110 V，供电电流1.5 A. 图 3为试验振荡器产生的逆M序列信号波形、功率谱和循环自相关时间序列.

图  3  逆M序列信号波形、功率谱和循环自相关时间序列

Figure  3.  Signal waveform, power spectrum and circular autocorrelation of inverse repeated M sequence

下载: 全尺寸图片 幻灯片

由于本项工作是试验性的，在台站未敷设新的观测线路，只有借用工作的线路在正常观测的间隙进行试验，且由于试验样机的参数是通过模拟得到，需在实际观测中经不断地调整和改正来逐步完善，所以观测进行了一天. 2016年7月30日9时开始试验，由于整点台站的其它仪器工作，所以从每个整点后30分钟开始工作，记录时长为20分钟，试验工作至17时50分结束；用1000 Hz采样频率同时记录激发电流信号和接收电压信号.从记录的电压曲线(图 4)看出，测量信号峰峰电压可达3500 mV(各种干扰及开关瞬时电磁感应、介质响应等综合值)，观测的编码信号地电介质响应电压为500 mV (图 4中黑点)，这两个信号的强度相差7倍，所以利用传统的地电观测方法，当地下介质的电性发生微弱变化时，基本上被噪声信号所淹没，想检测出目标层电阻率的真实变化则比较困难，而利用编码源的信号互相关理论，可以有效压制不相关噪声，分辨出有用信息，分辨信号能力较直流电法有一定的改善.

图  4  兰州观象台记录的测量电压(空心圆)和加载的编码源信号的电压(黑点)

Figure  4.  The value of ambient noise (open circles) and the encoding pseudorandom signals (dark dots) recorded at Lanzhou observatory

下载: 全尺寸图片 幻灯片

2.2.2   地电阻率计算

对图 4的电压信号按图 1的流程使用上述数据处理方法，产生参考信号，然后通过循环互相关法计算输入、输出信号与参考信号的互相关函数，再转换至频率域获得阻抗频率响应，乘以装置系数K即可得到不同频率的地电阻率，结果如图 5所示.可见：在不同时段、同一地点测量的地电阻率变化形态相同；以频率1 Hz为分界线，频率小于1 Hz的地电阻率变化相对稳定，地电阻率在6—10 Ω·m范围内变化，且随着频率的减小，地电阻率也在变小，反映了地下介质的电性结构在深部是低阻层，这与该台站电性分层一致(李乐进，1985)：表层为第四纪覆盖层(亚砂土、粉砂、沙砾石)，随着深度增加下伏为砂岩. 20世纪70年代该地区无电磁干扰测量电测曲线为Q型(钱复业等，1987)，根据该台站电测深曲线反演获得上部3层地电断面的参数：h₁=3.5 m, ρ₁= 80 Ω·m；h₂= 15 m, ρ₂=30 Ω·m；h₃=300 m，ρ₃=9 Ω·m.随频率增大，基于1 Hz以上频率点计算的地电阻率数值变大，分析认为有两个原因：其一，随着频率的增大，反映表层的信息更加丰富；其二，供电线路与接收线路在一起平行布线产生了电磁耦合. 图 6给出了2016年7月30日兰州观象台在不同频率点(0.05—1 Hz)对应时段记录的地电阻率曲线，可以看出低频点(0.05—0.2 Hz)的地电阻率与直流观测数值基本接近(图 2)，且受电磁耦合影响小，反映出地电阻率随时间的变化，且随着频率的变化，电阻率的变化层次也十分清晰，勾勒出不同频率电阻率的分布和变化，为研究目标层的电性精细变化提供了基础.

图  5  兰州观象台在不同时段记录0.01—25 Hz频率点的地电阻率曲线

Figure  5.  Georesistivity spectrum during 0.01-25 Hz at different time observed at Lanzhou observatory

下载: 全尺寸图片 幻灯片

图  6  2016年7月30日兰州观象台在不同频率点(0.01—1.03 Hz)对应时段记录的地电阻率曲线

Figure  6.  The georesistivity histogram at different time moment recorded at Lanzhou obser-vatory on July 30, 2016 for different frequencies (0.01-1.03 Hz)

下载: 全尺寸图片 幻灯片

2.3   陇南汉王台试验

陇南汉王地电台地处南北地震带中段，武都山字型弧状顶端转折部位.测区地形北高南低, 西高东低，北靠山、南临白龙江，覆盖厚度由东至西、由南向北从约22 m增至100 m以上，表层为碎砾石沙土，下伏中上志留纪碳硅质板岩、含黄铁矿片岩和磁黄铁矿砂岩互层，测区位于磁异常区附近，电测深曲线为Q型^①. 1974年建台布设3个地电阻率测道: N85°E道和N54°W道(AB=1 km，MN=0.333 km), N73°W道(AB=0.924 km，MN=0.308 km).近年来该地电测区由于受到铁路和高速公路、电站及工程建设的影响，观测资料质量大大降低. 图 7为2013—2014年N85°E测道的地电阻整点值曲线. 2013年7月开始由于电站的工作附加高频干扰，8月底A供电极附近填土建设电阻率曲线上升一个台阶；由于兰渝铁路汉王段通过汉王地电台的测区，其中铁路路基横跨N85°E道和N73°W道，共用供电电极(B₁和B₃)需要搬迁，故在2013年将该电极移至铁路的同侧路基下方；2014年2月初开始，铁路建设开始筑土碾压路基，在地面加高了3—4 m，在此过程中电阻率出现了持续振荡变化；由于测区环境严重干扰，2015年地面观测被迫停测.本次试验就选用2016年10月3—4日N85°E测道旧线路在原场地进行，以检验仪器的抗干扰性和数据观测质量.

① 国家地震局兰州地震研究所台群处编. 1983.甘肃省地震台志(中).

图  7  2013—2014年陇南汉王台N85°E测道地电阻率观测曲线

Figure  7.  The georesistivity curve recorded by the channel of N85°E at Hanwang observatory in 2013-2014

下载: 全尺寸图片 幻灯片

2.3.1   输入信号和输出电压

为了更利于压制干扰，通过仿真优选，使产生的信号序列周期是工频的整数倍，来设置编码源信号产生参数：寄存器位数m=6，序列周期126位；晶振频率6.3 MHz，分频后获得15 Hz时钟频率激励移位寄存器，产生的逆重复M序列周期为16.8 s.直流电压110 V，供电电流1.8 A.

2016年10月3日9时开始试验，从每个整点后20分开始工作，记录时间30分钟，工作至晚上8时结束；使用2个采集器按1000 Hz采样频率记录输入电流和输出电压信号.

2.3.2   地电阻率计算

地电阻率计算按照式(1)—(6) 计算.从2016年10月3日和4日两天的地电阻率观测(图 8)可看出，除10月3日10时的地电阻率测值变化较大外，其它时段测值比较稳定，10月3日0.06 Hz地电阻率在21.3 Ω·m上下波动，较平稳，10月4日0.06 Hz地电阻率在21.35 Ω·m上下波动，更为平稳.其它频率显示出相似的特征，高频点由于电磁耦合阻值升高，1.07 Hz的地电阻率在24.68 Ω·m上下波动.对10月3—4日0.06—1.07 Hz频率范围内电阻率变化的分析结果(表 1)表明，10月3日地电阻率的均值、变化幅度、变化率(日变化幅值与日观测均值的比值乘以100%)均大于4日，其原因在于3日数据采集没有同步，4日同步采样后变幅明显减小，各频率点的日变化率均小于1%.随着频率减小，地电阻率也减小，反映了陇南汉王地电测区地层的电性变化，表层为碎砾石沙土覆盖层，下伏含黄铁矿片岩和磁黄铁矿砂岩互层碳硅质板岩低阻层，与20世纪70年代测区无电测干扰测量电测深曲线Q型变化形态相吻合.

图  8  2016年10月3日(a)和4日(b)陇南汉王台在不同时段记录不同频率点的地电阻率曲线

Figure  8.  The georesistivity curves at different frequency at Hanwang observatory on October 3 (a) and 4 (b), 2016

下载: 全尺寸图片 幻灯片

表  1  2016年10月3日和4日陇南汉王台测试地电阻率变化幅值及变化率

Table  1.  The amplitude and variation of georesistivity ratio on October 3 and 4, 2016 at Hanwang observatory, Longnan city

频率/Hz	2016年10月3日				2016年10月4日
	ρs/(Ω·m)	Δρs/(Ω·m)	Δρs /ρs		ρs/(Ω·m)	Δρs/(Ω·m)	Δρs /ρs
0.06	21.43	0.65	3.03		21.38	0.15	0.70
0.19	21.50	0.25	1.16		21.47	0.21	0.98
0.31	21.69	0.42	1.92		21.69	0.09	0.41
0.44	21.98	0.36	1.63		22.02	0.17	0.77
0.57	22.40	0.55	2.45		22.42	0.15	0.67
0.69	22.79	1.09	4.76		22.91	0.12	0.52
0.82	23.32	1.26	5.39		23.44	0.19	0.81
0.94	23.86	1.56	6.54		24.06	0.21	0.87
1.07	24.44	2.13	8.72		24.72	0.18	0.73

下载: 导出CSV 

| 显示表格

3.   讨论与结论

甘肃省兰州观象台和陇南汉王地震台站的试验验证了本文所述的编码源电阻率测量原理是可行的，利用编码源地电仪在原有线路和场地可观测到高质量的地电阻率数据.尽管试验时间短，不能刻画出地电阻率的长期变化，但是通过试验得到的数据说明因干扰停测的汉王地电观测台，用编码源地电仪可以在原测区、原线路下工作，这说明改进的观测仪器方案是能够推广的.试验工作取得的成果如下：

1) 将相关辨识技术应用于地电阻率观测，可极大地压制观测场地的电磁干扰，提高地电阻率观测精度；

2) 由于采用编码源宽频信号激发，可得到一系列频率的地电阻率谱，利用这些地电阻率谱以期获得其它反映地震孕育电性变化的指标参数；

3) 通过两个地电台站试验，形成了地电阻率高精度观测方案，由辨识结果实时处理得到地电阻率随时间的变化形态，且随着频率的变化勾勒出不同频率电阻率的分布，对微小的变化具有较强的分辨力, 为研究目标层的电性变化提供了基础；

4) 获得了一些优化参数，运用这些参数改进编码源仪器可观测到高质量数据.两次试验结果表明，编码源产生参数如序列周期长度、叠加周期数、采样率等对观测结果均有明显影响.若序列周期是工频信号周期的偶数倍，将有助于提升抗干扰能力，故选序列长度以6阶寄存器产生为宜，时钟频率使产生的序列周期须在9—36 s范围内选取.

后续工作将通过仪器研发和观测试验完善观测系统.由于待监测地电系统非常复杂，地质情况也复杂，观测装置及布极的深度位置和极距大小均会影响到电流场在地电系统中的扩散路径，因此利用地震台站的现有观测装置，进一步研究不同深度电极装置对辨识结果的影响，对于地电阻率高精度观测有现实意义.

进一步研究去电磁耦合问题.电磁耦合是地电观测普遍存在的问题，电磁耦合去耦好坏会对地电阻率的观测精度产生直接的影响，利用编码源宽域广谱的特点，一次供电可获得一系列频率的地电阻率谱，为电磁去耦提供了条件(可用多频去耦，该方法后续撰文介绍)；相位信息也能较灵敏地反映地电阻率变化，所以本次试验对相位的观测进行了尝试，需要收发装置的采样频率严格同步，才能满足观测要求.目前通过方法技术创新研发的编码源地电阻率观测仪器已在装配，所形成的技术方法可以为利用现有观测装置进行地震台站地电阻率观测提供技术支持和实践基础.

本研究在实施过程中得到了中国地震局地质研究所赵国泽研究员、预测研究所钱家栋研究员、北京大学黄清华教授、长江大学胡文宝教授的关注, 得到了甘肃省地震局兰州观象台和陇南地震台的大力协助, 两位匿名审稿专家对本文提出了建设性的修改意见和建议，作者在此一并表示感谢.