CN102426393A - 电法勘探方法和装置 - Google Patents

Abstract

本发明涉及一种电法勘探方法，其通过发送机发送信号以及接收机对接收到的信号进行处理，得到冲激响应和频率响应，其特征在于，所述方法包括以下步骤：所述发送机发送伪随机序列，并利用所述伪随机序列对所述接收机接收到的时间序列进行互相关运算，进而得到时域冲激响应；将所述时域冲激响应变换到频域，从而得到频率响应；利用双科尔-科尔模型进行拟合和反演，得到复电阻率法的参数。

Description

电法勘探方法和装置

技术领域

本发明涉及电法勘探领域，具体而言，涉及一种电法勘探方法和装置。

背景技术

在寻找金属矿产方面电法勘探(特别是激电法)的应用非常广泛，但是由于生产矿山及其附近存在大量的高压输电线和各种用电设施，这将产生极其强烈的电磁干扰信号；地下开采坑道内大量的金属管道、车辆轨道和废石渣都给电法的供电和电场接收造成很大的困难。所有这一切，都影响了电法勘探方法的应用效果，因此电法勘探在解决危机矿山问题时存在严重的噪声干扰问题。

生产矿山附近电磁干扰尽管很大，但一般是随机性的干扰，而相关辨识技术是一种可以有效地去除随机噪声干扰的系统辨识方法，在机械、自动化、仪器仪表、认知科学、生物信息学等领域都已经得到了广泛的应用，但是在地球物理探测方面却没有得到足够的重视和应用。

目前，电法勘探中一般是通过堆叠、滤波、远参考等方法去除噪声，但是在矿山电法勘探中这些方法用于去除噪声就显得无能为力，因此勘察效果差，直接带来人力财力的巨大浪费。而且，为了提高信噪比，一个常用的方法就是提高发送电源的功率，这必然造成发送机异常笨重，这给山区的勘探工作带来巨大的困难。

另一方面，激电法是电法勘探中用于找矿的最有效的手段，特别是复电阻率法，能够提供较多的电性参数，探测精度较高。复电阻率法是指通过测量得到复电阻率谱也就是频率响应的一种电法勘探方法，但是需要测量多个频点的电场响应才能得到一条频率响应曲线，而且低频的时候需要的时间很长，因此探测效率很低，这是造成复电阻率法效果虽好却不能得到广泛运用的原因。

在电法勘探中如果采用相关辨识技术，那么一次测量就能得到一条频率响应曲线，可以有效解决复电阻率法的效率低下问题，并且具有复电阻率法探测精度高的优点，因此这是一种新的时间域测量的复电阻率法。

图1(a)和(b)示出了现有技术中复电阻率法的观测装置和等效电路。

传统的复电阻率法是通过在相当宽的超低频段上观测视复电阻率的幅度谱和相位谱或实分量和虚分量谱，以研究地下介质情况。如图1(a)和(b)所示，发送电流源AB每次发送一个频率的正弦波或者方波，接收端子MN测得对应的电压信号，多次测量得到一个频段内的多个频率的幅度和相位特性，得到复电阻率谱，再利用科尔-科尔模型进行拟合和反演，得到直流电阻率、充电率、时间常数和频率相关系数等电性参数进行地球物理解释。

这种方法能提供比较丰富的激电信息，但需在许多频率上做观测才能获得较完整的频谱，所以生产效率低。

发明内容

本发明的目的是利用相关辨识技术去除随机噪声；在时间域测量一次就可以得到一条频率响应曲线从而提高探测效率；另外，通过采用双科尔-科尔模型进行解释，不仅可以去除电磁耦合效应的影响，还可以得到电磁参数，探测精度提高。

根据本发明的一方面，提供一种电法勘探方法，其利用测量装置进行勘探，所述测量装置包括发送机和接收机，其特征在于，所述方法包括以下步骤：所述发送机发送伪随机序列v(t)，接收机接收时间序列信号u(t)；利用所述伪随机序列v(t)对所述接收机接收到的时间序列u(t)进行互相关运算，得到时域冲激响应h_e(t)；将所述时域冲激响应h_e(t)变换到频域，从而得到频率响应H_e(ω)；利用双科尔-科尔模型进行拟合和反演，得到复电阻率法的4个激电参数ρ₀、m₁、c₁、τ₁以及电磁耦合效应的3个参数m₂、c₂、τ₂；其中，ρ₀表示零频视电阻率：反映电极排列勘探体积内的平均电阻率；m₁为激电效应的视充电率：激电效应强度参数(％)，与电极排列勘探体积内的可极化物质体积含量正相关；τ₁为激电效应的视时间常数：激电效应特征参数(秒)，与电极排列勘探体积内的可极化物质的粒度大小等结构信息相关；c₁为激电效应的视频率相关系数，是无量纲的基地啊效应过程参数；与电极排列勘探体积内的可极化物质的激电效应类型以及极化物质混合分布均匀性相关；m₂为电磁效应的视充电率；τ₂为电磁效应的视时间常数；c₂为电磁效应的视频率相关系数。

优选地，所述方法还包括步骤：执行反卷积运算去除所述测量装置的影响。

优选地，通过多极距的观测方式得到目标极化体的真谱参数和几何分布。

优选地，所述伪随机序列为m序列或逆重复m序列。

优选地，将电极、导线、接地阻抗及测量装置的影响归结为一个系统，它的冲激响应为h_s(t)，大地系统的冲激响应记为h_e(t)，则

u(t)＝v(t)*h_e(t)*h_s(t)+n(t)

两边用输入伪随机序列v(t)进行互相关运算，消去随机噪声的影响，得到

R_vu(t)＝R_vv(t)*h_e(t)*h_s(t)

因此，先求输入和输出信号的互相关R_vu(t)以及输入信号的自相关R_vv(t)，反卷积运算后得到观测系统和大地系统总的冲激响应h_e(t)*h_s(t)，然后再通过反卷积运算去除观测系统h_s(t)的影响得到大地系统的冲激响应h_e(t)；

或者在频域实现上述运算：

P_vu(ω)＝P_vv(ω)·H_e(ω)·H_s(ω)

得到待辨识系统的频率响应： $H_e\left(\mathrm{\ω}\right)=\frac{P_{\mathrm{vu}}\left(\mathrm{\ω}\right)}{P_{\mathrm{vv}}\left(\mathrm{\ω}\right)\·H_s\left(\mathrm{\ω}\right)}$

将在发送机附近相关辨识得到的结果就认为是观测系统的冲激响应h_s(t)，变到频域即为H_s(ω)；

计算复电阻率谱ρ(ω)：

ρ(ω)＝K·H_e(ω)

其中K称为装置系数，

$K=\frac{2\mathrm{\π}}{\frac{1}{\mathrm{AM}}-\frac{1}{\mathrm{AN}}+\frac{1}{\mathrm{BN}}-\frac{1}{\mathrm{BM}}}$

式中AM、AN、BN、BM为电极间的距离，其中A、B是一对发送电极，M、N是一对接收电极；

利用双科尔-科尔模型

$\mathrm{\ρ}\left(\mathrm{\ω}\right)={\mathrm{\ρ}}_0\{1-m_1[1-\frac{1}{1+{\left(j{\mathrm{\ω\τ}}_1\right)}^{c_1}}\left]\right\}\{1-m_2[1-\frac{1}{1+{\left(j{\mathrm{\ω\τ}}_2\right)}^{c_2}}\left]\right\}$

拟合得到4个激电参数ρ₀、m₁、c₁、τ₁，以及3个电磁耦合效应参数m₂、c₂、τ₂。

优选地，拟合时ρ₀、m₁、c₁、τ₁和m₂、c₂、τ₂初值的选择方法如下：c₁的取值范围在0.1-0.6之间，c₂的取值范围为0.9-1.0，m₁、m₂的取值范围在0-0.98之间，而且m₁＜m₂，ρ₀的初值设为幅度谱或者实分量谱的低频渐近线的值，根据虚分量谱在高频段的峰值频率求出τ₂的初值，τ₁的初值定为τ₂的初值的50倍。

根据本发明的另一方面，提供一种用于电法勘探的发送机，其特征在于所述发送机包括全球定位系统GPS(Global Positioning System)模块(1)、伪随机信号发生器(2)、驱动及保护模块(3)、智能功率模块IPM(Intelligent Power Modules)(5)、温度控制模块(4)和供电模块(7)，和发送电极AB(6)；其中，GPS模块(1)通过串口向伪随机信号发生器(2)传递授时同步信息，完成数据信息的存储，并控制驱动及保护模块(3)；驱动及保护模块(3)驱动IPM(5)工作，保护电路用于保护整个系统的安全；温度控制系统(4)监测发射机内部的温度，当超过预定值后复位机器或关闭机器；IPM(5)将输出信号转换成用于电法勘探的电流信号并向发送电极AB(6)输送电流；供电模块(7)为板卡提供低压电源，并为IPM(5)提供所需高压。

优选地，伪随机信号发生器(2)进一步根据需要输出不同阶数、不同频率的m序列或者逆重复m序列，配置LCD用于显示发射机信号类型、参数及状态。

根据本发明的另一方面，提供一种用于电法勘探的接收机，其特征在于所述接收机包括GPS模块(8)、伪随机信号发生器(9)、接收电极MN(10)、前置放大模块(11)、滤波及陷波模块(12)、精密放大模块(13)、自相关器(14)、互相关器(15)、第一反卷积器(16)、第二反卷积器(17)、FFT变换器(18)、复电阻率谱生成器(19)、拟合器(20)、反演器(21)；其中，GPS模块(8)通过串口向伪随机信号发生器(9)传递授时同步信息，实现精准同步并且和发送机的伪随机序列保持精确同相，伪随机信号发生器(9)根据需要输出与发送机相同阶数、相同频率的m序列；接收电极MN(10)接收通过地下介质传来的电位差，然后经过前置放大模块(11)将接收到的电磁信号进行初步的放大，并经过滤波及陷波模块(12)滤除高频噪声及系统的工频噪声，再通过精密放大模块(13)将滤除噪声的信号进行低噪精密放大；自相关器(14)实现伪随机信号的自相关运算；互相关器(15)实现伪随机信号和接收信号的互相关运算；经过第一反卷积器(16)得到总的冲激响应，第二反卷积器(17)去除观测系统影响处理后获得大地系统冲激响应，FFT变换器(18)对大地冲激响应进行处理可获得大地系统的频率响应，大地系统频率响应经过复电阻率谱生成器(19)得到大地系统的复电阻率谱，拟合器(20)依据双科尔-科尔模型拟合可以获得视谱参数，反演器(21)进行联合反演获得真谱参数，从而获得地球物理有用的参数信息，并进行实时数据图形显示及存储。

根据本发明的再一方面，提供一种电法勘探系统，其特征在于包含上述发送机以及接收机。

与现有技术相比，根据本发明的技术方案能够带来以下的技术效果：

1、能够去除接收信号中叠加的随机干扰。根据维纳-霍夫方程，用输入的伪随机信号对输出信号进行互相关运算，由于输入信号和输出信号中叠加的随机噪声不相关，所以随机噪声被消除。

2、能够去除观测系统对测量结果的影响。对于一个实际系统，电极、导线、接地阻抗及测量装置都会对测量结果产生影响而降低辨识的准确度，利用反卷积运算，这个影响很容易被去除。

3、能够得到目标极化体的真谱参数和几何分布。单输入单输出辨识得到的冲激响应和频率响应反映的是发送接收电极之间勘探体积内地下介质的整体效应，采用多极距的高密度几何测深方法，能够得到目标极化体的真谱参数和几何分布，具有水平和纵向分辨力高、反映异常幅度大等优点，特别是在复杂地形条件下地质勘查效果较好。

4、能够去除电磁耦合效应的影响。激电法中如何去除电磁耦合效应的影响一直以来都是人们亟待解决的问题。利用双科尔-科尔模型的复电阻率法不仅可以得到激电效应的4个参数，还可以得到关于电磁耦合效应的3个参数，因此不仅可以分离电磁效应对激电参数的影响，利用电磁参数还可以得到更多的信息，因此探测精度大大提高。

附图说明

图1(a)示出了复电阻率法的观测装置；

图1(b)示出了复电阻率法的等效电路；

图2示出了根据本发明的电法勘探方法的辨识模型；

图3示出了待辨识系统的冲激响应(以一个典型的二阶系统为例)；

图4(a)示出了一个周期的m序列；

图4(b)示出了一个周期的逆重复m序列；

图4(c)示出了高斯白噪声的示意图；

图5示出了利用m序列辨识得到的冲激响应和真实冲激响应对比示意图；

图6示出了利用逆重复m序列辨识得到的冲激响应和真实冲激响应对比示意图；

图7示出了根据本发明的电法勘探方法的辨识模型中考虑观测系统的影响的情形；

图8示出了待辨识系统冲激响应反卷积前后的对比示意图；

图9示出了多极距观测系统的示意图；

图10示出了双科尔-科尔模型幅度谱辨识前后的对比示意图；

图11示出了双科尔-科尔模型相位谱辨识前后的对比示意图；

图12示出了双科尔-科尔模型拟合的幅度/相位示意图；

图13示出了双科尔-科尔模型拟合的实部/虚部示意图；

图14示出了根据本发明的电法勘探装置中发送机的示意图；

图15示出了根据本发明的电法勘探装置中接收机的示意图。

具体实施方式

图2示出了根据本发明的电法勘探方法的辨识模型。它的基本原理是维纳-霍夫方程。它的核心是采用伪随机序列(包括m序列和逆重复m序列)对接收到的时间序列进行相关处理，得到冲激响应，把冲激响应变到频域就是频率响应。

系统输入为v(t)，输出u(t)也就是测量到的信号，h(t)是待辨识系统的冲激响应，r(t)是输入v(t)经过被测系统后的响应，外部的随机噪声干扰n(t)，其中输入v(t)、输出u(t)是已知的，干扰n(t)为不可测的、未知的。

系统输入输出信号v(t)、u(t)有如下关系：

u(t)＝r(t)+n(t)＝v(t)*h(t)+n(t)    (1)

式中h(t)为被测系统的时域冲激响应，*为卷积运算。

式(1)两边同时用输入信号v(t)进行互相关运算：

R_vu(t)＝R_vv(t)*h(t)+R_vn(t)         (2)

式中R_vu(t)表示输入和输出的互相关，R_vv(t)表示输入的自相关，R_vn(t)表示干扰n(t)和输入信号的互相关。

由于外部的随机噪声干扰和输入信号不相关，因此(2)式中的R_vn(t)为0。可以看到，互相关运算的结果是消去了随机噪声的干扰。因此得到

R_vu(t)＝R_vv(t)*h(t)                (3)

这就是著名的维纳-霍夫方程。因此，我们只要求出输入的自相关函数R_vv(t)和输入输出的互相关R_vu(t)，就可以通过解卷积的方法得到被测系统的时域冲激响应h(t)。变换到频域就是频率响应。这样一次发送和接收就得到了一条频率响应曲线，包括幅频响应和相频响应，大大提高了探测效率。

但是如果输入是一般信号，那么求解方程(3)是比较困难的。为此我们要寻找一些特殊的信号作为输入信号以简化求解过程。而白噪声信号的自相关函数具有特殊的形式，其自相关函数是一个δ函数。当被辨识系统的输入为白噪声时，只要求出输入输出信号之间的互相关函数，就可以相除直接得到系统的冲激响应。但白噪声只是数学上的一个抽象，工程上不容易实现，更不具有可重复性。工程中常使用伪随机信号作为输入信号。它具有近似白噪声的性质，所以可以保证良好的辨识精度，而且在工程上容易实现。

辨识通常采用的伪随机信号是m序列或者逆重复m序列。m序列是指最大长度移位寄存器序列。逆重复m序列也称L序列，是对m序列扩展一倍，然后隔位取反后产生的，由于在一个周期内均值为0，因此在进行系统相关辨识时具有比m序列更优良的特性。附图3是一个待辨识系统的冲激响应h(t)，以一个典型的二阶系统为例。附图4(a)是一个周期，长度为31的m序列，图4(b)是一个周期的逆重复m序列，图4(c)是叠加的噪声n(t)。附图5是利用m序列进行相关辨识的结果，可以看到，经过相关辨识恢复的冲激响应基本和真实的冲激响应重合，噪声的影响变得很小。附图6是利用逆重复m序列进行相关辨识的结果，辨识效果比m序列更好。

另外，在本发明的一个优选实施例中，为了去除电极、导线、接地阻抗及测量装置的影响，在接收机的数据处理过程中进行了反卷积运算，其原理如下：

将电极、导线、接地阻抗及测量装置的影响归结为一个系统，它的冲激响应为h_s(t)，大地系统的冲激响应记为h_e(t)，如附图7所示。那么(1)式变成

u(t)＝v(t)*h_e(t)*h_s(t)+n(t)    (4)

(4)式两边用输入伪随机序列进行互相关运算，就可以消去随机噪声的影响，得到

R_vu(t)＝R_vv(t)*h_e(t)*h_s(t)     (5)

因此，先求输入和输出信号的互相关R_vu(t)以及输入信号的自相关R_vv(t)，反卷积运算后得到观测系统和大地系统总的冲激响应h_e(t)*h_s(t)，然后再通过反卷积运算去除观测系统h_s(t)的影响得到大地系统的冲激响应h_e(t)。从附图8可以看出，去除观测系统的影响后的冲激响应更接近待辨识系统的冲激响应。

也可以把卷积和反卷积运算转化成在频域或者复频域相乘或者相除的运算。因此(5)式变为：

P_vu(ω)＝P_vv(ω)·H_e(ω)·H_s(ω)    (6)

这样就得到待辨识系统的频率响应：

$H_e\left(\mathrm{\ω}\right)=\frac{P_{\mathrm{vu}}\left(\mathrm{\ω}\right)}{P_{\mathrm{vv}}\left(\mathrm{\ω}\right)\·H_s\left(\mathrm{\ω}\right)}---\left(7\right)$

在发送机附近相关辨识得到的结果就认为是观测系统的冲激响应h_s(t)，因为在发送机附近，可以忽略大地系统的影响，认为只有观测系统的影响。把它变到频域即为H_s(ω)。

在地球物理中用到的是复电阻率谱ρ(ω)。它和频率响应的关系是：

ρ(ω)＝K·H_e(ω)    (8)

这里K称为装置系数。如图1(a)所示，A、B是一对发送电极，M、N是一对接收电极，那么装置系数的大小取决于A、M、N、B四电极的相对位置：

$K=\frac{2\mathrm{\π}}{\frac{1}{\mathrm{AM}}-\frac{1}{\mathrm{AN}}+\frac{1}{\mathrm{BN}}-\frac{1}{\mathrm{BM}}}---\left(9\right)$

式中AM、AN、BN、BM为电极间的距离。

得到复电阻率谱以后，利用复电阻率法来进行地球物理解释。通常复电阻率法采用单科尔-科尔模型进行拟合：

$\mathrm{\ρ}\left(\mathrm{\ω}\right)={\mathrm{\ρ}}_0\{1-m_1[1-\frac{1}{1+{\left(j{\mathrm{\ω\τ}}_1\right)}^{c_1}}\left]\right\}---\left(10\right)$

这样拟合得到4个激电参数ρ₀、m、τ、c。由于复电阻率谱中含有电磁耦合效应的影响，所以这样得到的激电参数是不够准确的。

而利用双科尔-科尔模型(式(11)所示)进行拟合不仅可以得到4个激电参数ρ₀、m₁、c₁、τ₁，还可以得到关于电磁耦合效应的3个参数m₂、c₂、τ₂。

$\mathrm{\ρ}\left(\mathrm{\ω}\right)={\mathrm{\ρ}}_0\{1-m_1[1-\frac{1}{1+{\left(j{\mathrm{\ω\τ}}_1\right)}^{c_1}}\left]\right\}\{1-m_2[1-\frac{1}{1+{\left(j{\mathrm{\ω\τ}}_2\right)}^{c_2}}\left]\right\}---\left(11\right)$

这7个参数反映的是电极排列勘探体积内的总体信息，称为视谱参数，它们的物理意义是：

①ρ₀-表示零频视电阻率：反映电极排列勘探体积内的平均电阻率。

②m₁-激电效应的视充电率：激电效应强度参数(％)，与电极排列勘探体积内的可极化物质体积含量正相关。

③τ₁-激电效应的视时间常数：激电效应特征参数(秒)，与电极排列勘探体积内的可极化物质的粒度大小等结构信息相关。

④c₁-激电效应的视频率相关系数：激电效应过程参数(无量纲)，与电极排列勘探体积内的可极化物质的激电效应类型以及极化物质混合分布均匀性相关。

⑤m₂-电磁效应的视充电率。

⑥τ₂-电磁效应的视时间常数。

⑦c₂-电磁效应的视频率相关系数。

因此该方法不仅可以分离激电谱和电磁耦合谱，利用得到的电磁参数还能得到剩余电磁效应(REM)参数和电磁视电阻率ρ_ω等，因此可以提高探测精度。

在一个测点得到的是视谱参数，我们需要的是目标极化体(也就是需探明的矿体或者油气藏等)的真科尔-科尔参数和它的几何分布(包括深度、厚度、产状等信息)，所以需要2个以上的测点的信息。采取多极距观测，信息量就极大丰富了。因此在本发明的一个优选实施例中，采用了几何测深的方法——多极距观测，可以得到不同深度的异常体的信息。

这种工作方式如图9所示，AB是发送电极，M₁ N₁是第一对接收电极，N₁M₃(即为M₂ N₂)是第二对接收电极，M₃ N₃是第三对接收电极。同线排列的多个接收器同时接收发送电偶极源的电场响应，收发距按一定的间距增加。由于不同的收发距探测的深度不同，收发距越大反映的深度越大，因此在一个供电点一次供电，就可以接收到不同深度的多个频率的电场响应。沿剖面按一定的点距多次移动观测系统，就可观测到一条剖面上不同收发距的多个频率的径向电场响应，获取丰富的电场信息。可以采取以下两种方法由多极距观测结果得到目标极化体的真谱参数和几何分布。

①对多极距的测量结果也就是多个测点的视谱进行联合反演，可以求出目标极化体的几何分布及真谱参数。先对每个测点的复电阻率谱(包括振幅和相位)用常规复电阻率法反演分离出极化体和围岩的视谱，再对分离出的多个不同极距的视相位谱联合反演得到目标极化体的真谱参数和几何分布。

②利用分离出的电磁效应的参数可以计算剩余电磁效应(REM)参数和电磁视电阻率2个电磁参数，它们比常规视电阻率能更灵敏地反应地下导电性异常。这样反演每一个测点的复电阻率谱可获得六个参数：四个激电参数，两个电磁参数。多极距观测结果用六个参数的拟断面图表示，可以反映地电构造沿剖面和随深度的变化。

本领域的技术人员应当理解，在图1、图9中，发送端不一定用电偶极源，也可以使用磁性源线圈等任意一种电法勘探用主动源。接收端也不一定用偶极子接收，可以使用磁棒等任何一种电法勘探用接收器。另外，发送接收不一定采用如图9的多极距观测装置，可以是这种同线排列的，也可以是轴向测量，可以是这种一个发多个收的，也可以一个发一个收，或者多发多收。另外，对得到的频率响应，除了用复电阻率法来解释，还可以用直流电阻率法、激电法、电磁法来解释。在频率接近0时，频率响应曲线在该点的值就是该测点的直流视电阻率，得到多个测点的值后可以用直流电阻率法反演。在高频(例如大于10⁴Hz)，可以求取一些电磁参数，采用电磁法的资料处理和反演方法。在低频段则可以得到一些激电参数，例如视频散率、视频散率道间变化率、道间低频复电阻率的变化率等参数。另外，也可以直接对时间域的冲激响应进行地球物理解释。在考虑复电阻率法时，也不一定是这种乘积形式的双科尔-科尔模型，也可以是和形式的双科尔-科尔模型，还可以是Dias模型、Brown模型或者其他模型。

与现有技术相比较，根据本发明的技术方案的不同之处至少在于：

1、采用伪随机序列(包括m序列和逆重复m序列)对电法接收机接收到的时间序列进行互相关运算，得到系统的冲激响应和频率响应；

2、采用反卷积去除观测系统的影响；

3、采用多极距的观测方式得到目标极化体的真谱参数和几何分布；

4、对得到的频率响应采用复电阻率法进行反演解释。采用双科尔-科尔模型可以得到电磁参数，探测更精细。为了对本发明技术方案的技术效果有更好的了解，以下通过一个仿真实例进行说明。

由于典型二阶系统的特性是我们所熟知的，辨识效果也容易比对，所以下面以对一个二阶系统的辨识结果为例。假设一个待辨识二阶系统冲激响应为h(t)＝a·e^bt·sin(ct)，取a＝11.547，b＝-5，c＝8.66，如图3所示。该二阶系统截止频率为2Hz，它的冲激响应在1s后基本趋于零，所以系统的过渡时间约为1s。根据参数选择公式(12)和(13)，伪随机序列的码片宽度可以选择0.02s，而序列的周期长度为63即可。

$0.443\frac{1}{\mathrm{\Δ}}\≥f_{\max }---\left(12\right)$

$N=(1.2-1.5)\frac{T_5}{\mathrm{\Δ}}---\left(13\right)$

这里Δ为码片宽度，f_max为系统截止频率，T_s为系统的过渡时间，N为序列的周期长度。为了得到更好的辨识效果，取码片宽度为0.01s，序列长度取100个点。实际辨识过程中，可以发送几个周期(例如3-5个周期)的伪随机序列以得到更好的辨识效果。接收机也可以选择效果较好的波形进行运算。

在对实际的地下介质系统辨识的过程中，为了兼顾高频和低频分辨度以及和移位寄存器个数的矛盾，可以分成几个频段进行。另外，为了提高辨识精度可以提高采样频率，在一个码片宽度内可以采集多个点。

高斯白噪声是常见的噪声类型，图4(c)是采样频率为1kHz，10000个采样点的高斯白噪声，噪声的标准差为1。图4(a)是发送电极发出的伪随机m序列电流信号，图4(b)是发送电极发出的伪随机逆重复m序列电流信号，图5和图6分别是用m序列和逆重复m序列相关辨识的结果。可以看到，经过相关辨识以后，随机噪声对辨识得到的冲激响应影响很小，而且，采用逆重复m序列辨识效果更好。

考虑去除观测系统的影响。图8是待辨识的典型二阶系统的实际冲激响应、反卷积前后的冲激响应的对比，可以看到，利用反卷积去除观测系统的影响以后的结果(白色的实线)更接近实际系统的冲激响应(红色的+字线)，两者基本重合，而没有经过反卷积的结果(绿色的菱形曲线)则与实际系统相差较大。

这样得到的冲激响应经过傅里叶变换变到频域，就是H_e(ω)。利用(8)式得到复电阻率谱ρ(ω)。

对于一个双科尔-科尔模型系统(以乘积形式为例)，假定设定的参数分别为ρ₀＝25，m₁＝0.1，c₁＝0.12，τ₁＝1，m₂＝0.67，c₂＝0.95，τ₂＝0.01，加入了标准差为1V的高斯白噪声，相关辨识结果复电阻率谱ρ(ω)可以表示成幅度谱(见图10)和相位谱(见图11)。辨识后的结果(白色+字线)只取了正频率部分。可以看到辨识效果较好，基本和实际系统谱重合。

采用经典的反演拟合方法阻尼最小二乘法对得到的复电阻率谱进行拟合，可以得到ρ₀，m₁，c₁，τ₁，m₂，c₂，τ₂共7个激电和电磁参数，参见式(11)。

采用阻尼最小二乘法进行拟合时初值的选取是非常关键的，初值选择不当可能会造成不收敛的情况，因而得不到需要的结果。对双科尔-科尔模型进行拟合共有7个待定参数，参数多拟合更不容易得到理想的效果，因此选择合适的初值就显得更加关键。这里详细说明初值的选择方法。

首先，这几个参数是有取值范围的：c₁的取值范围在0.1-0.6之间，c₂的取值范围为0.9-1.0，m₁、m₂的取值范围在0-0.98之间，而且m₁＜m₂，τ₁＞＞τ₂。而零频电阻率的值ρ₀就是幅度或者实部的低频渐近线的值。

τ₂也很容易得到。在高频段电磁耦合占绝对优势的时候，虚分量有个峰值，这个峰值频率就对应电磁耦合效应的时间常数τ₂。再根据τ₁＞＞τ₂可以估计τ₁。

m₁、c₁和m₂、c₂的变化范围都不大，因此在它们的取值范围内取个估计值就可以了，比如m₁＝0.25、c₁＝0.25、m₂＝0.75、c₂＝1。

另外，为了防止迭代过程中因为修改步长太大导致最后的迭代结果不在参数的有效范围内，甚至出现负值，在程序中可以预先设定ρ₀、m₁、c₁、τ₁、m₂、c₂、τ₂七个参数的取值范围，既可以加快迭代拟合的速度又可以确保拟合的准确性。

图12和13是对一个双科尔-科尔模型(参数设定为ρ₀＝25，m₁＝0.1，c₁＝0.12，τ₁＝0.1，m₂＝0.67，c₂＝0.95，τ₂＝0.003)的拟合结果，迭代8次便达到了设定的误差门限要求。实践证明，这种方法是行之有效的。

上述过程是对一个测点得到该测点的视谱参数的情况。为了得到目标极化体的真谱参数和几何分布，需要进行多点观测。多极距观测方式是一种高密度的几何测深方式，这种工作方式如图9所示，在一个供电点一次供电，就可以接收到不同深度的不同频率的电场响应。AB是发送电极，M₁ N₁是第一对接收电极，N₁M₃(即为M₂ N₂)是第二对接收电极，M₃ N₃是第三对接收电极。利用(8)式得到不同测点也就是不同极距下的复电阻率谱，再利用阻尼最小二乘法反演拟合得到7个视谱参数，然后再利用联合反演方法得到真谱参数和几何分布，或者是由这7个视谱参数得到6个激电、电磁参数的拟断面图，表示地电构造沿剖面和随深度的变化。

图14和图15分别示出了根据本发明的电法勘探装置中发送机和接收机的示意图。要注意的是，这两图仅是通过功能模块的方式示出了其中的主要模块，并不表示具体的物理结构模块。

在图14中，GPS模块(1)通过串口向伪随机信号发生器(2)传递授时同步信息，伪随机信号发生器(2)可根据需要输出不同阶数、不同频率的m序列或者逆重复m序列，配置LCD用于显示发射机信号类型、参数及状态，同时伪随机信号发生器(2)完成数据信息的存储，并控制驱动及保护模块(3)。驱动及保护模块(3)中驱动电路用于将产生的伪随机序列驱动IPM(5)工作，保护电路用于保护整个系统的安全。温度控制系统(4)时刻监测发射机内部的温度，当超过预定值后复位机器甚至关闭机器。IPM(5)将输出信号转换成用于电法勘探的电流信号并向发送电极AB(6)输送电流。供电模块(7)为板卡提供低压电源，并为IPM(5)提供所需高压。

在图15中，GPS模块(8)通过串口向伪随机信号发生器(9)传递授时同步信息，实现精准同步并且和发送机的伪随机序列保持精确同相，伪随机信号发生器(9)可根据需要输出与发送机相同阶数、相同频率的m序列。接收电极MN(10)接收通过地下介质传来的电位差，然后经过前置放大模块(11)将接收到的电磁信号进行初步的放大，并经过滤波及陷波模块(12)滤除高频噪声及系统的工频噪声，再通过精密放大模块(13)将滤除噪声的信号进行低噪精密放大。自相关器(14)实现伪随机信号的自相关运算。互相关器(15)实现伪随机信号和接收信号的互相关运算。自相关器(14)实现伪随机信号的自相关运算；互相关器(15)实现伪随机信号和接收信号的互相关运算；经过第一反卷积器(16)得到总的冲激响应，第二反卷积器(17)去除观测系统影响处理后获得大地系统冲激响应，FFT变换器(18)对大地冲激响应进行处理可获得大地系统的频率响应，大地系统频率响应经过复电阻率谱生成器(19)得到大地系统的复电阻率谱，拟合器(20)依据科尔-科尔模型拟合可以获得视谱参数，反演器(21)进行联合反演获得的真谱参数，从而获得地球物理有用的参数信息，并可以进行实时数据图形显示及存储。

Claims (10)

1.一种电法勘探方法，其利用测量装置进行勘探，所述测量装置包括发送机和接收机，其特征在于，所述方法包括以下步骤：

(1)所述发送机发送伪随机序列v(t)，接收机接收时间序列信号u(t)；

(2)利用所述伪随机序列v(t)对所述接收机接收到的时间序列u(t)进行互相关运算，得到时域冲激响应h_e(t)；

(3)将所述时域冲激响应h_e(t)变换到频域，从而得到频率响应H_s(ω)；

(4)利用双科尔-科尔模型进行拟合和反演，得到复电阻率法的4个激电参数ρ₀、m₁、c₁、τ₁以及电磁耦合效应的3个参数m₂、c₂、τ₂；

其中，ρ₀表示零频视电阻率：反映电极排列勘探体积内的平均电阻率；m₁为激电效应的视充电率：激电效应强度参数(％)，与电极排列勘探体积内的可极化物质体积含量正相关；τ₁为激电效应的视时间常数：激电效应特征参数(秒)，与电极排列勘探体积内的可极化物质的粒度大小等结构信息相关；c₁为激电效应的视频率相关系数，是无量纲的基地啊效应过程参数；与电极排列勘探体积内的可极化物质的激电效应类型以及极化物质混合分布均匀性相关；m₂为电磁效应的视充电率；τ₂为电磁效应的视时间常数；c₂为电磁效应的视频率相关系数。

2.根据权利要求1所述的电法勘探方法，其特征在于，所述方法还包括步骤：执行反卷积运算去除所述测量装置的影响。

3.根据权利要求1或2所述的电法勘探方法，其特征在于，通过多极距的观测方式得到目标极化体的真谱参数和几何分布。

4.根据权利要求1或3所述的电法勘探方法，其特征在于，所述伪随机序列为m序列或逆重复m序列。

5.根据权利要求1或4所述的电法勘探方法，其中将电极、导线、接地阻抗及测量装置的影响归结为一个系统，它的冲激响应为h_s(t)，大地系统的冲激响应记为h_e(t)，则

u(t)＝v(t)*h_e(t)*h_s(t)+n(t)

两边用输入伪随机序列v(t)进行互相关运算，消去随机噪声的影响，得到

R_vu(t)＝R_vv(t)*h_e(t)*h_s(t)

因此，先求输入和输出信号的互相关R_vu(t)以及输入信号的自相关R_vv(t)，反卷积运算后得到观测系统和大地系统总的冲激响应h_e(t)*h_s(t)，然后再通过反卷积运算去除观测系统h_s(t)的影响得到大地系统的冲激响应h_e(t)；

或者在频域实现上述运算：

P_vu(ω)＝P_vv(ω)·H_e(ω)·H_s(ω)

得到待辨识系统的频率响应：

$H_e\left(\mathrm{\ω}\right)=\frac{P_{\mathrm{vu}}\left(\mathrm{\ω}\right)}{P_{\mathrm{vv}}\left(\mathrm{\ω}\right)\·H_s\left(\mathrm{\ω}\right)}$

将在发送机附近相关辨识得到的结果就认为是观测系统的冲激响应h_s(t)，变到频域即为H_s(ω)；

计算复电阻率谱ρ(ω)；

ρ(ω)＝K·H_e(ω)

其中K称为装置系数，

$K=\frac{2\mathrm{\π}}{\frac{1}{\mathrm{AM}}-\frac{1}{\mathrm{AN}}+\frac{1}{\mathrm{BN}}-\frac{1}{\mathrm{BM}}}$

式中AM、AN、BN、BM为电极间的距离，其中A、B是一对发送电极，M、N是一对接收电极；

利用双科尔-科尔模型

$\mathrm{\ρ}\left(\mathrm{\ω}\right)={\mathrm{\ρ}}_0\{1-m_1[1-\frac{1}{1+{\left(j{\mathrm{\ω\τ}}_1\right)}^{c_1}}\left]\right\}\{1-m_2[1-\frac{1}{1+{\left(j{\mathrm{\ω\τ}}_2\right)}^{c_2}}\left]\right\}$

拟合得到4个激电参数ρ₀、m₁、c₁、τ₁，以及3个电磁耦合效应参数m₂、c₂、τ₂。

6.根据权利要求5所述的电法勘探方法，其中拟合时ρ₀、m₁、c₁、τ₁和m₂、c₂、τ₂初值的选择方法如下：c₁的取值范围在0.1-0.6之间，c₂的取值范围为0.9-1.0，m₁、m₂的取值范围在0-0.98之间，而且m₁＜m₂，ρ₀的初值设为幅度谱或者实分量谱的低频渐近线的值，根据虚分量谱在高频段的峰值频率求出τ₂的初值，τ₁的初值定为τ₂的初值的50倍。

7.一种用于电法勘探的发送机，其特征在于所述发送机包括全球定位系统GPS(Global Positioning System)模块(1)、伪随机信号发生器(2)、驱动及保护模块(3)、智能功率模块IPM(Intelligent Power Modules)(5)、温度控制模块(4)和供电模块(7)，和发送电极AB(6)；

其中，GPS模块(1)通过串口向伪随机信号发生器(2)传递授时同步信息，完成数据信息的存储，并控制驱动及保护模块(3)；驱动及保护模块(3)驱动IPM(5)工作，保护电路用于保护整个系统的安全；温度控制系统(4)监测发射机内部的温度，当超过预定值后复位机器或关闭机器；IPM(5)将输出信号转换成用于电法勘探的电流信号并向发送电极AB(6)输送电流；供电模块(7)为板卡提供低压电源，并为IPM(5)提供所需高压。

8.根据权利要求7所述的用于电法勘探的发送机，其特征在于伪随机信号发生器(2)进一步根据需要输出不同阶数、不同频率的m序列或者逆重复m序列，配置LCD用于显示发射机信号类型、参数及状态。

9.一种用于电法勘探的接收机，其特征在于所述接收机包括GPS模块(8)、伪随机信号发生器(9)、接收电极MN(10)、前置放大模块(11)、滤波及陷波模块(12)、精密放大模块(13)、自相关器(14)、互相关器(15)、第一反卷积器(16)、第二反卷积器(17)、FFT变换器(18)、复电阻率谱生成器(19)、拟合器(20)、反演器(21)；

其中，GPS模块(8)通过串口向伪随机信号发生器(9)传递授时同步信息，实现精准同步并且和发送机的伪随机序列保持精确同相，伪随机信号发生器(9)根据需要输出与发送机相同阶数、相同频率的伪随机序列；接收电极MN(10)接收通过地下介质传来的电位差，然后经过前置放大模块(11)将接收到的电磁信号进行初步的放大，并经过滤波及陷波模块(12)滤除高频噪声及系统的工频噪声，再通过精密放大模块(13)将滤除噪声的信号进行低噪精密放大；自相关器(14)实现伪随机信号的自相关运算；互相关器(15)实现伪随机信号和接收信号的互相关运算；经过第一反卷积器(16)得到总的冲激响应，第二反卷积器(17)去除观测系统影响处理后获得大地系统冲激响应，FFT变换器(18)对大地冲激响应进行处理获得大地系统的频率响应，大地系统频率响应经过复电阻率谱生成器(19)得到大地系统的复电阻率谱，拟合器(20)依据双科尔-科尔模型拟合获得视谱参数，反演器(21)进行联合反演获得真谱参数，从而获得地球物理参数信息，并进行实时数据图形显示及存储。

10.一种电法勘探系统，其特征在于包含如权利要求7或8所述的发送机以及如权利要求9所述的接收机。

Images