CN105607131A - 编码源电磁测深法获取大地全区视电阻率谱的方法及装置 - Google Patents

Abstract

本发明提供了一种编码源电磁测深法获取大地全区视电阻率谱的方法及装置。该方法包括：向待测大地输入电流信号，采集待测大地输出的电磁场响应信号时间序列，电流信号为逆重复M序列且周期为工频周期的偶数倍；根据电流信号和电磁场响应信号时间序列，通过循环互相关方式获取待测大地的频率响应；根据待测大地的频率响应及接收电磁场分量，获取待测大地的全区视电阻率谱。本发明通过逆重复M序列的电流信号激发，采用循环互相关技术获取全区视电阻率谱，提高了获取全区视电阻率的效率，及提高了获取的全区视电阻率的频率分辨率，且由于采用互相关法及电流信号逆重复M序列的周期为工频周期的偶数倍，能有效压制工频干扰，增强了抗干扰能力。

Description

编码源电磁测深法获取大地全区视电阻率谱的方法及装置

技术领域

本发明涉及信号采集与处理领域、系统辨识领域及地球物理电磁测深法勘探领域，具体而言，涉及一种编码源电磁测深法获取大地全区视电阻率谱的方法及装置。

背景技术

由于不同频率的电磁场在大地中的趋肤深度不同，因此在电磁测深领域中可以通过在地表测量电磁场分量，获取待测大地的全区视电阻率，根据全区视电阻率随频率的变化，推测待测大地的电性深度分布情况。

当前，相关技术中通常采用扫频法来获取待测大地的视电阻率。在应用扫频法时，预先设置多个不同的频率，在勘测时依次向待测大地输入预设的多个不同频率的电流信号，根据待测大地对多个不同频率的电流信号对应的电磁响应输出信号进行观测，获得待测大地的视电阻率。

由于依次向待测大地输入预设的多个不同频率的电流信号，如此一个频率一个频率的进行扫频勘测，效率很低，且没有采用相关辨识技术，抗干扰能力很低。

发明内容

有鉴于此，本发明实施例的目的在于提供一种编码源电磁测深法获取大地全区视电阻率谱的方法及装置，提高了获取全区视电阻率的效率，及提高了获取的全区视电阻率的频率分辨率，且能有效压制工频干扰，增强了抗干扰能力。

第一方面，本发明实施例提供了一种编码源电磁测深法获取大地全区视电阻率谱的方法，所述方法包括：

向待测大地输入电流信号，采集所述待测大地输出的电磁场响应信号时间序列，所述电流信号为逆重复M序列且周期为工频周期的偶数倍；

根据所述电流信号和所述电磁场响应信号时间序列，通过循环互相关方式获取所述待测大地的频率响应；

根据所述待测大地的频率响应及接收电磁场分量，获取所述待测大地的全区视电阻率谱。

结合第一方面，本发明实施例提供了上述第一方面的第一种可能的实现方式，其中，

所述根据所述电流信号和所述电磁场响应信号时间序列，通过循环互相关方式获取所述待测大地的频率响应，包括：

按照所述电流信号的时间序列周期和采样率触发软件程序生成参考信号，将所述电磁场响应信号时间序列及所述电流信号分别与所述参考信号进行循环互相关，得到所述电磁场响应信号时间序列对应的互相关时间序列及所述电流信号对应的互相关时间序列；

对所述电磁场响应信号时间序列对应的互相关时间序列及所述电流信号对应的互相关时间序列分别进行快速傅立叶变换，得到所述电磁场响应信号时间序列的互功率谱及所述电流信号的互功率谱；

根据所述电磁场响应信号时间序列的互功率谱及所述电流信号的互功率谱，获取所述待测大地的频率响应。

结合第一方面，本发明实施例提供了上述第一方面的第二种可能的实现方式，其中，

所述接收电磁场分量为水平电场时间序列；

所述根据所述待测大地的频率响应及接收电磁场分量，获取所述待测大地的全区视电阻率谱，包括：

根据所述待测大地的频率响应，通过公式(1)获得所述待测大地的阻抗谱与接收电磁场分量及大地电阻率的关系式；

对公式(1)对应的所述待测大地的阻抗谱对应的关系式进行变形，得到公式(2)所示的所述待测大地的全区视电阻率谱；

$Z_{Ex}\left(\mathrm{\ω}\right)=\frac{\mathrm{\ρ}\left(\mathrm{\ω}\right)dL\·MN}{2{\mathrm{\πr}}^3}{C}_{Ex}\mathrm{...}\left(1\right)$

${\mathrm{\ρ}}_{Ex}\left(\mathrm{\ω}\right)=\frac{2{\mathrm{\πr}}^3\·Z_{Ex}\left(\mathrm{\ω}\right)}{dL\·MN\·C_{Ex}}\mathrm{...}\left(2\right)$

其中，在公式(1)和(2)中，Z_Ex(ω)为大地阻抗，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，MN为接收电极的电极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离，为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角，i为虚数单位，为均匀大地的波数，μ₀为自由空间的磁导率，ρ_Ex(ω)为所述全区视电阻率。

结合第一方面，本发明实施例提供了上述第一方面的第三种可能的实现方式，其中，

所述接收电磁场分量为垂直磁场时间序列；

所述根据所述待测大地的频率响应及接收电磁场分量，获取所述待测大地的全区视电阻率谱，包括：

根据所述待测大地的频率响应，通过公式(3)获得所述待测大地的频率响应与接收电磁场分量及大地电阻率的关系式；

对公式(3)对应的所述频率响应对应的关系式进行变形，得到公式(4)所示的所述待测大地的全区视电阻率谱；

其中，在公式(3)和(4)中，Z_Hz(ω)为所述频率响应，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离，μ₀为自由空间的磁导率，i为虚数单位，为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角， $C_{Hz}=1-(1-{\mathrm{ik}}_{1}r-\frac{1}{3}{\left(k_{1}r\right)}^2)\exp \left({\mathrm{ik}}_{1}r\right),k_1^2={\mathrm{i\ω\μ}}_0/\mathrm{\ρ}$ 为均匀大地的波数，ρ_Hz(ω)为所述全区视电阻率。

结合第一方面，本发明实施例提供了上述第一方面的第四种可能的实现方式，其中，所述方法还包括：

根据所述待测大地的频率响应，获取所述待测大地的传递函数，由所述传递函数计算所述待测大地的阶跃响应或冲激响应，根据所述阶跃响应的晚期渐近值或所述冲激响应的峰值时刻，计算所述待测大地的视电阻率。

第二方面，本发明实施例提供了一种编码源电磁测深法获取大地全区视电阻率谱的装置，所述装置包括：

采集模块，用于向待测大地输入电流信号，采集所述待测大地输出的电磁场响应信号时间序列，所述电流信号为逆重复M序列且周期为工频周期的偶数倍；

第一获取模块，用于根据所述电流信号和所述电磁场响应信号时间序列，通过循环互相关方式获取所述待测大地的频率响应；

第二获取模块，用于根据所述待测大地的频率响应及接收电磁场分量，获取所述待测大地的全区视电阻率谱。

结合第二方面，本发明实施例提供了上述第二方面的第一种可能的实现方式，其中，

所述第一获取模块包括：

循环互相关单元，用于按照所述电流信号的时间序列周期和采样率触发软件程序生成参考信号，将所述电磁场响应信号时间序列及所述电流信号分别与所述参考信号进行循环互相关，得到所述电磁场响应信号时间序列对应的互相关时间序列及所述电流信号对应的互相关时间序列；

傅立叶变换单元，用于对所述电磁场响应信号时间序列对应的互相关时间序列及所述电流信号对应的互相关时间序列分别进行快速傅立叶变换，得到所述电磁场响应信号时间序列的互功率谱及所述电流信号的互功率谱；

获取单元，用于根据所述电磁场响应信号时间序列的互功率谱及所述电流信号的互功率谱，获取所述待测大地的频率响应。

结合第二方面，本发明实施例提供了上述第二方面的第二种可能的实现方式，其中，

所述接收电磁场分量为水平电场时间序列；

所述第二获取模块包括：

第一获得单元，用于根据所述待测大地的频率响应，通过公式(1)获得所述待测大地的阻抗谱与接收电磁场分量及大地电阻率的关系式；

第一计算单元，用于对公式(1)对应的所述待测大地的阻抗谱对应的关系式进行变形，得到公式(2)所示的所述待测大地的全区视电阻率谱；

$Z_{Ex}\left(\mathrm{\ω}\right)=\frac{\mathrm{\ρ}\left(\mathrm{\ω}\right)dL\·MN}{2{\mathrm{\πr}}^3}{C}_{Ex}\mathrm{...}\left(1\right)$

${\mathrm{\ρ}}_{Ex}\left(\mathrm{\ω}\right)=\frac{2{\mathrm{\πr}}^3\·Z_{Ex}\left(\mathrm{\ω}\right)}{dL\·MN\·C_{Ex}}\mathrm{...}\left(2\right)$

其中，在公式(1)和(2)中，Z_Ex(ω)为大地阻抗，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，MN为接收电极的电极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离，为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角，i为虚数单位，为均匀大地的波数，μ₀为自由空间的磁导率，ρ_Ex(ω)为所述全区视电阻率。

结合第二方面，本发明实施例提供了上述第二方面的第三种可能的实现方式，其中，

所述接收电磁场分量为垂直磁场时间序列；

所述第二获取模块包括：

第二获得单元，用于根据所述待测大地的频率响应，通过公式(3)获得所述待测大地的频率响应与接收电磁场分量及大地电阻率的关系式；

第二计算单元，用于对公式(3)对应的所述频率响应对应的关系式进行变形，得到公式(4)所示的所述待测大地的全区视电阻率谱；

其中，在公式(3)和(4)中，Z_Hz(ω)为所述频率响应，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离，μ₀为自由空间的磁导率，i为虚数单位，为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角， $C_{Hz}=1-(1-{\mathrm{ik}}_{1}r-\frac{1}{3}{\left(k_{1}r\right)}^2)\exp \left({\mathrm{ik}}_{1}r\right),k_1^2={\mathrm{i\ω\μ}}_0/\mathrm{\ρ}$ 为均匀大地的波数，ρ_Hz(ω)为所述全区视电阻率。

结合第二方面，本发明实施例提供了上述第二方面的第四种可能的实现方式，其中，所述装置还包括：

计算模块，用于根据所述待测大地的频率响应，获取所述待测大地的传递函数，由所述传递函数计算所述待测大地的阶跃响应或冲激响应，根据所述阶跃响应的晚期渐近值或所述冲激响应的峰值时刻，计算所述待测大地的视电阻率。

在本发明实施例提供的方法及装置中，向待测大地输入电流信号，采集待测大地输出的电磁场响应信号时间序列，电流信号为逆重复M序列且周期为工频周期的偶数倍；根据电流信号和电磁场响应信号时间序列，通过循环互相关方式获取待测大地的频率响应；根据待测大地的频率响应及接收电磁场分量，获取待测大地的全区视电阻率谱。本发明通过逆重复M序列的电流信号激发，采用循环互相关技术获取全区视电阻率谱，提高了获取全区视电阻率的效率，及提高了获取的全区视电阻率的频率分辨率，且由于采用互相关法及电流信号逆重复M序列的周期为工频周期的偶数倍，能有效压制工频干扰，增强了抗干扰能力。

为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1A示出了本发明实施例1所提供的一种获取全区视电阻率的方法流程图；

图1B示出了本发明实施例1所提供的一种逆重复M序列的时间域波形示意图；

图1C示出了本发明实施例1所提供的一种逆重复M序列的功率谱示意图；

图1D示出了本发明实施例1所提供的一种逆重复M序列的自相关函数曲线的示意图；

图1E示出了本发明实施例1所提供的一种编码源电磁测深阵列式工作布置示意图；

图1F示出了本发明实施例1所提供的一种编码源电磁测深发射系统的结构示意图；

图1G示出了本发明实施例1所提供的一种编码源电磁测深接收机系统的结构示意图；

图2示出了本发明实施例2所提供的一种获取全区视电阻率的装置结构示意图。

具体实施方式

下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

考虑到相关技术中依次向待测大地输入预设的多个不同频率的电流信号，如此一个频率一个频率的进行扫频勘测，效率很低，且没有采用相关辨识技术，抗干扰能力不强。基于此，本发明实施例提供了一种编码源电磁测深法获取大地全区视电阻率谱的方法及装置。下面通过实施例进行描述。

实施例1

参见图1A，本发明实施例提供了一种编码源电磁测深法获取大地全区视电阻率谱的方法。该方法具体包括以下步骤：

步骤101：向待测大地输入电流信号，采集待测大地输出的电磁场响应信号时间序列，电流信号为逆重复M序列且周期为工频周期的偶数倍。

本发明实施例中，通过时钟信号触发线性反馈移位寄存器生成电流信号，该电流信号为周期性伪随机的逆重复M序列。

本发明实施例的执行主体为用于编码源电磁测深法获取大地全区视电阻率谱的装置，该装置中设置有N+1位寄存器，N为大于或等于3的整数。在这N+1位寄存器中，1位寄存器用于生成二分频移位寄存器时钟信号，剩余的N位寄存器组成线性反馈移位寄存器，预先设定线性反馈移位寄存器的反馈位，这N位线性反馈移位寄存器组成了最大长度伪随机m序列产生电路。上述另1位寄存器二分频线性反馈移位寄存器时钟信号，将该时钟信号与N位线性反馈移位寄存器生成的伪随机m序列进行异或运算，实现对伪随机m序列隔位取反，得到逆重复M序列。驱动大功率IGBT(InsulatedGateBipolarTransistor，绝缘栅双极型晶体管)逆变器，通过水平电偶极子向大地供入编码电流。用接收电磁场相同的接收机，按照设定的采样率将输入大地的编码电流信号采样为时间序列，将采样得到的时间序列信号确定为电流信号，电流信号为多周期的信号。

其中，逆重复M序列的周期及频带是可调节的，可以通过选择不同的线性反馈移位寄存器位数及时钟频率来调节。线性反馈移位寄存器的位数一般取5、6、7、8、9或10位，时钟频率可以在频率区间[0.5Hz，10kHz]内进行选择。其中，如表1所示，不同的线性反馈移位寄存器位数与时钟频率，对应的逆重复M序列的周期及序列位长不同。表1中本原多项式为产生逆重复M序列的反馈位表达式，时钟频率与逆重复M序列周期的对应表中只是示意性地给出时钟频率，具体操作中时钟频率可以在频率区间[0.5Hz，10kHz]内进行选择。改变产生逆重复M序列的时钟频率，可方便地实现分频段激发，一般2-3次激发就能获得0.1Hz至10000Hz的频率范围的频率覆盖，提高了工作效率。

表1

其中，逆重复M序列的周期是最大线性反馈伪随机m序列的周期的2倍。取逆重复M序列的周期为工频周期的偶数倍，即电流信号的周期为工频周期的偶数倍。在工业用电较多的场景下，电流信号会受工业用电的影响。而电流信号的相关峰值表现为周期性的正负交替出现，因此取电流信号的周期为工频周期的偶数倍，可以有效的压制工频干扰，削弱工频干扰对电流信号的影响，进而提高获取的频率响应函数的准确性。

其中，逆重复M序列的时间域波形如图1B所示，逆重复M序列1个周期内波形是随机的。逆重复M序列的功率谱如图1C所示，其功率谱为离散的线状谱。逆重复M序列的自相关函数曲线如图1D所示，该自相关函数的曲线表现为正负交替的尖脉冲。

如图1E所示的编码源电磁测深阵列式工作布置方式中，通过水平电偶极子源AB向待测大地输入逆重复M序列的电流信号，并通过接收电极MN或磁棒采集待测大地输出的电磁场响应信号时间序列。在本发明实施例中，既可以同步方式采集电磁场响应信号时间序列，也可以异步方式采集电磁场响应信号时间序列。

步骤102：根据电流信号和电磁场响应信号时间序列，通过循环互相关方式获取待测大地的频率响应。

按照电流信号的时间序列周期和采样率触发软件程序生成参考信号，将电磁场响应信号时间序列及电流信号分别与参考信号进行循环互相关，得到电磁场响应信号时间序列对应的互相关时间序列及电流信号对应的互相关时间序列。

其中，该参考信号序列长度与参与计算的电磁场响应信号及输入电流信号序列长度相同，是输入电流信号周期的整数倍。该参考信号的生成方式可以与电流信号的生成方式相同，也可从电流信号的时间序列中提取。其中，设置产生逆重复M序列的线性反馈移位寄存器的位数及反馈位之后，产生的逆重复M序列是固定的，只是不同的寄存器初始状态产生的逆重复M序列之间具有时间延迟，但是产生的逆重复M序列的频谱分布是相同的。

其中，该参考信号的周期与电流信号周期相同，均为工频周期的偶数倍。后续根据该参考信号来对电流信号和电磁场响应信号进行循环互相关处理，可以有效压制工频干扰，提高抗干扰能力。

在本发明实施例中，根据上述参考信号和电流信号，通过如下公式(5)将参考信号与电流信号进行循环互相关，得到电流信号对应的互相关时间序列。以及，根据参考信号和电磁场响应信号时间序列，通过如下公式(6)将参考信号与电磁场响应信号时间序列进行循环互相关，得到电磁场响应信号时间序列对应的互相关时间序列。

C_us(t)＝cxcorr(u,SS)…(5)

C_ys(t)＝cxcorr(yout,SS)…(6)

其中，在公式(5)和(6)中，u为电流信号，yout为电磁场响应信号时间序列，SS为参考信号，t为时间，C_us(t)为电流信号对应的互相关时间序列，C_ys(t)为电磁场响应信号时间序列对应的互相关时间序列，cxcorr(u,SS)表示对电流信号u和参考信号SS进行循环互相关运算，cxcorr(yout,SS)表示对电磁场响应信号时间序列yout和参考信号SS进行循环互相关运算。

通过上述方式得到电磁场响应信号时间序列对应的互相关时间序列及电流信号对应的互相关时间序列后，对电磁场响应信号时间序列对应的互相关时间序列及电流信号对应的互相关时间序列分别进行快速傅立叶变换，得到电磁场响应信号时间序列的互功率谱及电流信号的互功率谱。

具体地，通过如下公式(7)对电流信号对应的互相关时间序列进行快速傅立叶变换，得到电流信号的互功率谱。以及通过如下公式(8)对电磁场响应信号时间序列对应的互相关时间序列进行快速傅立叶变换，得到电磁场响应信号时间序列的互功率谱。

$P_{C_{us}}\left(\mathrm{\ω}\right)=FFT\left(L_{us}\right(t\left)\right)\mathrm{...}\left(7\right)$

$P_{C_{ys}}\left(\mathrm{\ω}\right)=FFT\left(L_{ys}\right(t\left)\right)\mathrm{...}\left(8\right)$

其中，在公式(7)和(8)中，ω为频率，为电流信号的互功率谱，为电磁场响应信号时间序列的互功率谱，FFT(C_us(t))表示对电流信号对应的互相关时间序列C_us(t)进行快速傅立叶变换，FFT(C_ys(t))表示对电磁场响应信号时间序列对应的互相关时间序列C_ys(t)进行快速傅立叶变换。

通过上述方式得到电磁场响应信号时间序列的互功率谱及电流信号的互功率谱后，计算相同频率的电磁场响应信号时间序列的互功率谱与电流信号的互功率谱之间的比值，得到如公式(9)所示的待测大地的频率响应。

$H\left(\mathrm{\ω}\right)=\frac{P_{C_{ys}}\left(\mathrm{\ω}\right)}{P_{C_{us}}\left(\mathrm{\ω}\right)}\mathrm{...}\left(9\right)$

其中，在公式(9)中，ω为频率值，为频率值ω对应的电磁场响应信号时间序列的谱线峰值，为频率值ω对应的电流信号的谱线峰值，H(ω)为待测大地的频率响应。

在本发明实施例，当采用同步方式采集电流信号及电磁场响应信号时间序列的时间序列时，还根据电磁场响应信号时间序列的互功率谱及电流信号的互功率谱，获取待测大地的相位谱，待测大地的相位谱如公式(10)所示。

其中，在公式(10)中，mod为模运算符，为频率ω对应的电磁场响应信号时间序列的互功率谱的相位，为频率ω对应的电流信号的互功率谱的相位，为待测大地对应的相位谱。

另外，若采用异步方式采集电流信号及电磁场响应信号时间序列的时间序列时，只能获取待测大地的幅度谱，无法获得相位谱。

步骤103：根据待测大地的频率响应和接收电磁场分量，获取待测大地的全区视电阻率谱。

在本发明实施例中，接收电磁场分量可以为水平电场时间序列或垂直磁场时间序列，水平电场时间序列用符号Ex表示，垂直磁场时间序列用符号Hz表示。当接收电磁场分量不同时，获取待测大地的全区视电阻率谱的方式也不同。

当接收电磁场分量为水平电场时间序列时，通过如下方式获取待测大地的全区视电阻率谱，具体包括：

根据待测大地的频率响应，通过公式(1)获得所述待测大地的阻抗谱与接收电磁场分量及大地电阻率的关系式；对公式(1)对应的所述待测大地的阻抗谱对应的关系式进行变形，得到公式(2)所示的所述待测大地的全区视电阻率谱；

$Z_{Ex}\left(\mathrm{\ω}\right)=\frac{\mathrm{\ρ}\left(\mathrm{\ω}\right)dL\·MN}{2{\mathrm{\πr}}^3}{C}_{Ex}\mathrm{...}\left(1\right)$

${\mathrm{\ρ}}_{Ex}\left(\mathrm{\ω}\right)=\frac{2{\mathrm{\πr}}^3\·Z_{Ex}\left(\mathrm{\ω}\right)}{dL\·MN\·C_{Ex}}\mathrm{...}\left(2\right)$

其中，在公式(1)和(2)中，Z_Ex(ω)为大地阻抗，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，MN为接收电极的电极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离，为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角，i为虚数单位，为均匀大地的波数，μ₀为自由空间的磁导率，ρ_Ex(ω)为全区视电阻率。

在上述公式(1)和(2)中，dL、MN、r及均为观测装置的装置参数，根据场源及观测点位坐标确定这些参数，在观测到水平电场时间序列后，根据待测大地的频率响应、水平电场时间序列及这些装置参数，通过公式(1)和(2)来获得待测大地的全区视电阻率谱。

在本发明实施例中，当接收电磁场分量是水平电场时间序列Ex时，观测值是接收电极间的电压值时间序列，该电压值时间序列V_Ex＝E_x·MN，MN为接收电极的电极间距，所以由以上计算过程，公式(9)所示的频率响应可变形为如下公式(11)：

$Z_{Ex}\left(\mathrm{\ω}\right)=\frac{P_{{\mathrm{CV}}_{ex}s}\left(\mathrm{\ω}\right)}{P_{{\mathrm{CI}}_{AB}s}\left(\mathrm{\ω}\right)}\mathrm{...}\left(11\right)$

其中，在公式(11)中，为参考信号SS与电压值时间序列V_Ex的互功率谱，为参考信号SS与水平电偶源的两级发射的电流信号I_AB(t)的互功率谱。

由均匀大地表面水平电偶极子源产生的电场Ex的计算公式(12)：

其中，在公式(12)中ρ为均匀大地的电阻率，I是发射电流，r是接收电极的中点与水平电偶源的两极中点之间的距离，是接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角，dL是水平电偶源的两极间距，为均匀大地的波数，i为虚数单位，μ₀为自由空间的磁导率。

将上述公式(12)两边除以发射电流I，然后与公式(11)对应相等，即可得到公式(1)所示的待测大地的大地阻抗Z_Ex(ω)，其中Z_Ex(ω)是由观测电场Ex电压值和激发电流I循环互相关辨识得到的大地阻抗。通过对上述公式(1)进行变形即可得到如公式(2)所示的待测大地的全区视电阻率谱，其中公式(2)所示的全区视电阻率需要采用迭代法求解。

在本发明实施例中，当接收电磁场分量为垂直磁场时间序列Hz时，通过如下方式获取待测大地的全区视电阻率谱，具体包括：

根据待测大地的频率响应，通过公式(3)获得所述待测大地的频率响应与接收电磁场分量及大地电阻率的关系式；对公式(3)对应的所述频率响应对应的关系式进行变形，得到公式(4)所示的所述待测大地的全区视电阻率谱；

其中，在公式(3)和(4)中，Z_Hz(ω)为频率响应，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离，μ₀为自由空间的磁导率，i为虚数单位，为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角， $C_{Hz}=1-(1-{\mathrm{ik}}_{1}r-\frac{1}{3}{\left(k_{1}r\right)}^2)\exp \left({\mathrm{ik}}_{1}r\right),k_1^2={\mathrm{i\ω\μ}}_0/\mathrm{\ρ}$ 为均匀大地的波数，ρ_Hz(ω)为全区视电阻率。

在上述公式(3)和(4)中，dL、r及均为观测装置的装置参数，根据场源及观测点位坐标确定这些参数，在观测到垂直磁场时间序列后，根据待测大地的频率响应、垂直磁场时间序列及这些装置参数，通过公式(3)和(4)来获得待测大地的全区视电阻率谱。

在本发明实施例中，当接收电磁场分量是垂直磁场时间序列Hz时，以由以上计算过程，公式(9)所示的频率响应可变形为如下公式(13)：

$Z_{Hz}\left(\mathrm{\ω}\right)=\frac{P_{CHzs}\left(\mathrm{\ω}\right)}{P_{{\mathrm{CI}}_{AB}s}\left(\mathrm{\ω}\right)}\mathrm{...}\left(13\right)$

其中，在公式(13)中，P_CHzs(ω)为参考信号SS与垂直磁场分量时间序列Hz的互功率谱，为参考信号SS与水平电偶源的发射电流I_AB(t)的互功率谱。

由均匀大地表面水平电偶极子源产生的垂直磁场Hz计算公式(14)：

其中，在公式(14)中I是发射电流，r是磁场接收点至水平电偶源的两极中点之间的距离，是磁场接收点至水平电偶源的两极中点之间的连线与水平电偶源的两极延长线的夹角。dL是水平电偶源的两极间距。为均匀大地的波数，μ₀是自由空间的磁导率，i为虚数单位。

将上述公式(14)两边除以发射电流I，然后与公式(13)对应相等，即可得到公式(3)所示的待测大地的频率响应Z_Hz(ω)，其中Z_Hz(ω)是由观测垂直磁场Hz和激发电流I循环互相关辨识得到的大地频率域响应。通过对上述公式(3)进行变形即可得到如公式(4)所示的待测大地的全区视电阻率谱，其中公式(4)所示的全区视电阻率需要采用迭代法求解。

在本发明实施例中，除通过上述方式获取待测大地的全区视电阻率谱外，还可以通过辨识待测大地的传递函数来获取大地视电阻率，具体包括：

根据待测大地的频率响应，获取待测大地的传递函数，由传递函数计算待测大地的阶跃响应或冲激响应，根据阶跃响应的晚期渐近值或冲激响应的峰值时刻，计算待测大地的视电阻率。

因为不同频率的电磁场在大地中的趋肤深度不一样，所获得的全区视电阻率谱反映不同深度的电阻率变化，从而达到电磁测深的目的。

在本发明实施例中，由于所采用的激发信号逆重复M序列编码源是宽频带信号，所以通过本发明实施例提供的方法同时可得到一系列频率的视电阻率谱。为了提高获得的全区视电阻率谱的数据质量，需要发射和接收多个周期的逆重复M序列编码源信号，如此可以提高获取全区视电阻率的信噪比。

在本发明实施例中，通过如图1F所示的编码源电磁测深发射系统来发射电流信号，编码源电磁测深发射系统能选择输出编码源激励波形，能选择激发频带。激励波形包括寄存器位数不同的逆重复M序列，产生逆重复M序列的移位寄存器位数宜为5，6，7，8，9，10，11。移位寄存器位数越多，产生的逆重复M序列周期越长，信号频带越宽。产生的激发频带由激励移位寄存器动作的时钟频率确定，时钟频率越高，信号频带越高。

通过如图1G所示的编码源电磁测深接收机系统来接收电磁场响应信号时间序列。接收机系统能采集单道或多道电场信号和磁场信号。对于不同频带的电磁场信号采用不同的采样率采集时间序列，采样率应有多组可调，对于高频带信号采用较高的采样率采样，对于低频带信号采用较低的采样率采样。接收机系统配有能贮存大量数据的存贮空间及数据传输功能。编码源电磁测深接收机系统和编码源电磁测深发射系统可采用GPS(GlobalPositioningSystem，全球定位系统)同步或高精度石英钟同步。

在本发明实施例中，电流信号时间序列通过在发送供电回路中串接大功率取样电阻取样，用与接收电磁场相同的接收机，按相同的采样率采集激发电流时间序列。在本发明实施例中可采用轴向装置阵列式观测电场分量，或者采用赤道装置阵列式观测电场分量或垂直磁场分量均可高精度、高分辨率获得大地全区视电阻率谱。

在阵列式观测时，各接收机及接收电场或磁场传感器的一致性可通过仪器标定的方式解决。设标定的接收机频率响应为H_s(ω)，则由以下公式(15)可去除接收机系统自身的影响:

$H\left(\mathrm{\ω}\right)=\left(\frac{P_{Cys}\left(\mathrm{\ω}\right)}{P_{Cus}\left(\mathrm{\ω}\right)}\right)/H_s\left(\mathrm{\ω}\right)\mathrm{...}\left(15\right)$

为了提高观测精度和质量，采集系统配置质量监控系统，数据实时处理，在多周期电流激发下，观测精度达到设定要求方停止采样。

在本发明实施例中，向待测大地输入电流信号，采集待测大地输出的电磁场响应信号时间序列，电流信号为逆重复M序列且周期为工频周期的偶数倍；根据电流信号和电磁场响应信号时间序列，通过循环互相关方式获取待测大地的频率响应；根据待测大地的频率响应及接收电磁场分量，获取待测大地的全区视电阻率谱。本发明通过逆重复M序列的电流信号激发，采用循环互相关技术获取全区视电阻率谱，提高了获取全区视电阻率的效率，及提高了获取的全区视电阻率的频率分辨率，且由于采用互相关法及电流信号逆重复M序列的周期为工频周期的偶数倍，能有效压制工频干扰，增强了抗干扰能力。

实施例2

参见图2，本发明实施例提供了一种编码源电磁测深法获取大地全区视电阻率谱的装置，该装置用于执行上述实施例1提供的编码源电磁测深法获取大地全区视电阻率谱的方法。该装置具体包括：

采集模块201，用于向待测大地输入电流信号，采集待测大地输出的电磁场响应信号时间序列，电流信号为逆重复M序列且周期为工频周期的偶数倍；

第一获取模块202，用于根据电流信号和电磁场响应信号时间序列，通过循环互相关方式获取待测大地的频率响应；

第二获取模块203，用于根据待测大地的频率响应及接收电磁场分量，获取待测大地的全区视电阻率谱。

上述第一获取模块202通过循环互相关单元、傅立叶变换单元和获取单元来获取待测大地的频率响应。循环互相关单元，用于按照电流信号的时间序列周期和采样率触发软件程序生成参考信号，将电磁场响应信号时间序列及电流信号分别与参考信号进行循环互相关，得到电磁场响应信号时间序列对应的互相关时间序列及电流信号对应的互相关时间序列；傅立叶变换单元，用于对电磁场响应信号时间序列对应的互相关时间序列及电流信号对应的互相关时间序列分别进行快速傅立叶变换，得到电磁场响应信号时间序列的互功率谱及电流信号的互功率谱；获取单元，用于根据电磁场响应信号时间序列的互功率谱及电流信号的互功率谱，获取待测大地的频率响应。

当接收电磁场分量为水平电场时间序列时，第二获取模块203通过第一获得单元和第一计算单元来计算全区视电阻率。

第一获得单元，用于根据所述待测大地的频率响应，通过公式(1)获得所述待测大地的阻抗谱与接收电磁场分量及大地电阻率的关系式；

第一计算单元，用于对公式(1)对应的所述待测大地的阻抗谱对应的关系式进行变形，得到公式(2)所示的所述待测大地的全区视电阻率谱；

$Z_{Ex}\left(\mathrm{\ω}\right)=\frac{\mathrm{\ρ}\left(\mathrm{\ω}\right)dL\·MN}{2{\mathrm{\πr}}^3}{C}_{Ex}\mathrm{...}\left(1\right)$

${\mathrm{\ρ}}_{Ex}\left(\mathrm{\ω}\right)=\frac{2{\mathrm{\πr}}^3\·Z_{Ex}\left(\mathrm{\ω}\right)}{dL\·MN\·C_{Ex}}\mathrm{...}\left(2\right)$

其中，在公式(1)和(2)中，Z_Ex(ω)为大地阻抗，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，MN为接收电极的电极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离，为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角，i为虚数单位，为均匀大地的波数，μ₀为自由空间的磁导率，ρ_Ex(ω)为全区视电阻率。

在上述公式(1)和(2)中，dL、MN、r及均为观测装置的装置参数，根据场源及观测点位坐标确定这些参数，在观测到水平电场时间序列后，根据待测大地的频率响应、水平电场时间序列及这些装置参数，通过公式(1)和(2)来获得待测大地的全区视电阻率谱。

当接收电磁场分量为垂直磁场时间序列；

所述第二获取模块203包括：

第二获得单元，用于根据所述待测大地的频率响应，通过公式(3)获得所述待测大地的频率响应与接收电磁场分量及大地电阻率的关系式；

第二计算单元，用于对公式(3)对应的所述频率响应对应的关系式进行变形，得到公式(4)所示的所述待测大地的全区视电阻率谱；

其中，在公式(3)和(4)中，Z_Hz(ω)为频率响应，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离，μ₀为自由空间的磁导率，i为虚数单位，为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角， $C_{Hz}=1-(1-{\mathrm{ik}}_{1}r-\frac{1}{3}{\left(k_{1}r\right)}^2)\exp \left({\mathrm{ik}}_{1}r\right),k_1^2={\mathrm{i\ω\μ}}_0/\mathrm{\ρ}$ 为均匀大地的波数，ρ_Hz(ω)为全区视电阻率。

在上述公式(3)和(4)中，dL、r及均为观测装置的装置参数，根据场源及观测点位坐标确定这些参，在观测到垂直磁场时间序列后，根据待测大地的频率响应、垂直磁场时间序列及这些装置参数，通过公式(3)和(4)来获得待测大地的全区视电阻率谱。

在本发明实施例中，该装置还可以通过如下计算模块来获取全区视电阻率。

计算模块，用于根据所述待测大地的频率响应，获取所述待测大地的传递函数，由所述传递函数计算所述待测大地的阶跃响应或冲激响应，根据所述阶跃响应的晚期渐近值或所述冲激响应的峰值时刻，计算所述待测大地的视电阻率。

在本发明实施例中，向待测大地输入电流信号，采集待测大地输出的电磁场响应信号时间序列，电流信号为逆重复M序列且周期为工频周期的偶数倍；根据电流信号和电磁场响应信号时间序列，通过循环互相关方式获取待测大地的频率响应；根据待测大地的频率响应及接收电磁场分量，获取待测大地的全区视电阻率谱。本发明通过逆重复M序列的电流信号激发，采用循环互相关技术获取全区视电阻率谱，提高了获取全区视电阻率的效率，及提高了获取的全区视电阻率的频率分辨率，且由于采用互相关法及电流信号逆重复M序列的周期为工频周期的偶数倍，能有效压制工频干扰，增强了抗干扰能力。

本发明实施例所提供的编码源电磁测深法获取大地全区视电阻率谱的装置可以为设备上的特定硬件或者安装于设备上的软件或固件等。所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，前述描述的模块、装置和单元的具体工作过程，均可以参考上述方法实施例中的对应过程。

在本申请所提供的几个实施例中，应该理解到，所揭露装置和方法，可以通过其它的方式实现。以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，又例如，多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些通信接口，装置或单元的间接耦合或通信连接，可以是电性，机械或其它的形式。

所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。

另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。

所述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，RandomAccessMemory)、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (10)

1.一种编码源电磁测深法获取大地全区视电阻率谱的方法，其特征在于，所述方法包括：

向待测大地输入电流信号，采集所述待测大地输出的电磁场响应信号时间序列，所述电流信号为逆重复M序列且周期为工频周期的偶数倍；

根据所述电流信号和所述电磁场响应信号时间序列，通过循环互相关方式获取所述待测大地的频率响应；

根据所述待测大地的频率响应及接收电磁场分量，获取所述待测大地的全区视电阻率谱。

2.根据权利要求1所述的方法，其特征在于，所述根据所述电流信号和所述电磁场响应信号时间序列，通过循环互相关方式获取所述待测大地的频率响应，包括：

按照所述电流信号的时间序列周期和采样率触发软件程序生成参考信号，将所述电磁场响应信号时间序列及所述电流信号分别与所述参考信号进行循环互相关，得到所述电磁场响应信号时间序列对应的互相关时间序列及所述电流信号对应的互相关时间序列；

对所述电磁场响应信号时间序列对应的互相关时间序列及所述电流信号对应的互相关时间序列分别进行快速傅立叶变换，得到所述电磁场响应信号时间序列的互功率谱及所述电流信号的互功率谱；

根据所述电磁场响应信号时间序列的互功率谱及所述电流信号的互功率谱，获取所述待测大地的频率响应。

3.根据权利要求1所述的方法，其特征在于，所述接收电磁场分量为水平电场时间序列；

所述根据所述待测大地的频率响应及接收电磁场分量，获取所述待测大地的全区视电阻率谱，包括：

根据所述待测大地的频率响应，通过公式(1)获得所述待测大地的阻抗谱与接收电磁场分量及大地电阻率的关系式；

对公式(1)对应的所述待测大地的阻抗谱对应的关系式进行变形，得到公式(2)所示的所述待测大地的全区视电阻率谱；

$Z_{Ex}\left(\mathrm{\ω}\right)=\frac{\mathrm{\ρ}\left(\mathrm{\ω}\right)dL\·MN}{2{\mathrm{\πr}}^3}{C}_{Ex}...\left(1\right)$

${\mathrm{\ρ}}_{Ex}\left(\mathrm{\ω}\right)=\frac{2{\mathrm{\πr}}^3\·Z_{Ex}\left(\mathrm{\ω}\right)}{dL\·MN\·C_{Ex}}...\left(2\right)$

其中，在公式(1)和(2)中，Z_Ex(ω)为大地阻抗，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，MN为接收电极的电极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离， 为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角，i为虚数单位，为均匀大地的波数，μ₀为自由空间的磁导率，ρ_Ex(ω)为所述全区视电阻率。

4.根据权利要求1所述的方法，其特征在于，所述接收电磁场分量为垂直磁场时间序列；

所述根据所述待测大地的频率响应及接收电磁场分量，获取所述待测大地的全区视电阻率谱，包括：

根据所述待测大地的频率响应，通过公式(3)获得所述待测大地的频率响应与接收电磁场分量及大地电阻率的关系式；

对公式(3)对应的所述频率响应对应的关系式进行变形，得到公式(4)所示的所述待测大地的全区视电阻率谱；

其中，在公式(3)和(4)中，Z_Hz(ω)为所述频率响应，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离，μ₀为自由空间的磁导率，i为虚数单位，为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角， $C_{Hz}=1-(1-{\mathrm{ik}}_{1}r-\frac{1}{3}{\left(k_{1}r\right)}^2)\exp \left({\mathrm{ik}}_{1}r\right),k_1^2={\mathrm{i\ω\μ}}_0/\mathrm{\ρ}$ 为均匀大地的波数，ρ_Hz(ω)为所述全区视电阻率。

5.根据权利要求1所述的方法，其特征在于，所述方法还包括：

根据所述待测大地的频率响应，获取所述待测大地的传递函数，由所述传递函数计算所述待测大地的阶跃响应或冲激响应，根据所述阶跃响应的晚期渐近值或所述冲激响应的峰值时刻，计算所述待测大地的视电阻率。

6.一种编码源电磁测深法获取大地全区视电阻率谱的装置，其特征在于，所述装置包括：

采集模块，用于向待测大地输入电流信号，采集所述待测大地输出的电磁场响应信号时间序列，所述电流信号为逆重复M序列且周期为工频周期的偶数倍；

第一获取模块，用于根据所述电流信号和所述电磁场响应信号时间序列，通过循环互相关方式获取所述待测大地的频率响应；

第二获取模块，用于根据所述待测大地的频率响应及接收电磁场分量，获取所述待测大地的全区视电阻率谱。

7.根据权利要求6所述的装置，其特征在于，所述第一获取模块包括：

循环互相关单元，用于按照所述电流信号的时间序列周期和采样率触发软件程序生成参考信号，将所述电磁场响应信号时间序列及所述电流信号分别与所述参考信号进行循环互相关，得到所述电磁场响应信号时间序列对应的互相关时间序列及所述电流信号对应的互相关时间序列；

傅立叶变换单元，用于对所述电磁场响应信号时间序列对应的互相关时间序列及所述电流信号对应的互相关时间序列分别进行快速傅立叶变换，得到所述电磁场响应信号时间序列的互功率谱及所述电流信号的互功率谱；

获取单元，用于根据所述电磁场响应信号时间序列的互功率谱及所述电流信号的互功率谱，获取所述待测大地的频率响应。

8.根据权利要求6所述的装置，其特征在于，所述接收电磁场分量为水平电场时间序列；

所述第二获取模块包括：

第一获得单元，用于根据所述待测大地的频率响应，通过公式(1)获得所述待测大地的阻抗谱与接收电磁场分量及大地电阻率的关系式；

第一计算单元，用于对公式(1)对应的所述待测大地的阻抗谱对应的关系式进行变形，得到公式(2)所示的所述待测大地的全区视电阻率谱；

$Z_{Ex}\left(\mathrm{\ω}\right)=\frac{\mathrm{\ρ}\left(\mathrm{\ω}\right)dL\·MN}{2{\mathrm{\πr}}^3}{C}_{Ex}...\left(1\right)$

${\mathrm{\ρ}}_{Ex}\left(\mathrm{\ω}\right)=\frac{2{\mathrm{\πr}}^3\·Z_{Ex}\left(\mathrm{\ω}\right)}{dL\·MN\·C_{Ex}}...\left(2\right)$

其中，在公式(1)和(2)中，Z_Ex(ω)为大地阻抗，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，MN为接收电极的电极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离， 为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角，i为虚数单位，为均匀大地的波数，μ₀为自由空间的磁导率，ρ_Ex(ω)为所述全区视电阻率。

9.根据权利要求6所述的装置，其特征在于，所述接收电磁场分量为垂直磁场时间序列；

所述第二获取模块包括：

第二获得单元，用于根据所述待测大地的频率响应，通过公式(3)获得所述待测大地的频率响应与接收电磁场分量及大地电阻率的关系式；

第二计算单元，用于对公式(3)对应的所述频率响应对应的关系式进行变形，得到公式(4)所示的所述待测大地的全区视电阻率谱；

其中，在公式(3)和(4)中，Z_Hz(ω)为所述频率响应，ρ为均匀大地的电阻率，ω为频率，dL为水平电偶源的两极间距，r为接收电极的中点与水平电偶源的两极中点之间的距离，μ₀为自由空间的磁导率，i为虚数单位，为接收电极的中点至水平电偶源的两极中点之间的连线与水平电偶源的两极连线的夹角， $C_{Hz}=1-(1-{\mathrm{ik}}_{1}r-\frac{1}{3}{\left(k_{1}r\right)}^2)\exp \left({\mathrm{ik}}_{1}r\right),k_1^2={\mathrm{i\ω\μ}}_0/\mathrm{\ρ}$ 为均匀大地的波数，ρ_Hz(ω)为所述全区视电阻率。

10.根据权利要求6所述的装置，其特征在于，所述装置还包括：

计算模块，用于根据所述待测大地的频率响应，获取所述待测大地的传递函数，由所述传递函数计算所述待测大地的阶跃响应或冲激响应，根据所述阶跃响应的晚期渐近值或所述冲激响应的峰值时刻，计算所述待测大地的视电阻率。