CN105445805A - 一种时空阵列差分电磁勘探方法 - Google Patents

Abstract

本发明提供了一种时空阵列差分电磁勘探方法，包括以下步骤：(1)在空间多个测点及远参考点处进行多时窗的同步阵列数据采集；(2)根据观测数据构建测点信号时空数据矩阵X、天然电磁场信号时空数据矩阵X^r，以及相关噪声场信号时空数据矩阵Y；(3)对X^r进行矩阵分解求解天然电磁场源极化参数α，对Y进行矩阵分解求解相关噪声场源极化参数β，再利用α、β以及X求解测点对应于天然场源的空间模数U和对应于相关噪声场源的空间模数V；(4)利用U和V求解各个测点的天然场张量阻抗与相关噪声场张量阻抗。本发明数据采集方式灵活，且所有测道、时窗均允许存在噪声，根据场源信号的特征实现信噪分离，一次处理获得所有测点的天然场张量阻抗与相关噪声场张量阻抗。

Description

一种时空阵列差分电磁勘探方法

技术领域

本发明涉及一种勘察地球物理领域的电磁探测方法，特别涉及一种时空阵列差分电磁勘探方法。

背景技术

在勘察地球物理电磁法领域，天然场源电磁法(Magnetotelluric，MT；Audio-frequencyMagnetotelluric，AMT)的数据采集经历了标量测量、张量测量、连续剖面测量等发展阶段；尽管3D阵列测量目前已有研究与实践，但相应的数据处理技术仍然停留在单点处理阶段。相关噪声是制约天然场源电磁法在强干扰区应用的主要瓶颈之一，传统的最小二乘法、Robust估计、远参考法等均难以压制强相关噪声的影响。针对该问题，学者们先后发展了一系列处理手段，取得了一定的处理效果，但仍然存在相应的问题。如时间域数学形态滤波方法，通过提取含噪时间域波形进行信噪分离，时频域Hilbert-Huang变换处理方法，通过将信号在时频域进行分解并剔选进行信噪分离；这类方法在信噪识别、参数选取以及结果评价方面仍有待完善；基于远参考磁场相关性控制的功率谱数据删选方法以及基于远参考磁场与测点平面波磁场相互关系的信噪分离方法，通过远参考磁场提取测点数据中的平面波场信息，进而实现信噪分离；但这类方法对于强度大、持续时间长的相关噪声，以及远参考磁道中含有相关噪声的情况，处理效果仍然不能令人满意。

发明内容

本发明所解决的技术问题是，针对现有技术的不足，提供一种时空阵列差分电磁勘探方法，本发明所采用的阵列数据采集方式灵活，测点的空间阵列分布无需固定形式，基于统一的数据方程，尽可能的利用时空阵列中的所有观测数据，且所有测道、时窗均允许存在噪声，根据场源信号的特征实现信噪分离，一次处理获得所有同步阵列测点的天然场张量阻抗与相关噪声场张量阻抗，进而获得解释所需的地电参数。

本发明所采用的阵列数据采集方式灵活，测点的空间阵列分布无需固定形式，基于统一的数据方程，尽可能的利用时空阵列中的所有观测数据，且所有测道、时窗均允许存在噪声，根据场源信号的特征实现信噪分离，一次处理获得所有同步阵列测点的天然场张量阻抗与相关噪声场张量阻抗，进而获得解释所需的地电参数。

本发明的技术方案为：

一种时空阵列差分电磁勘探方法，包括以下步骤：

步骤1、观测设计：确定探测目标及勘探深度范围，在测区内设置测线及J个同步观测测点，J≥2；同时设置J'(J'≥1)个远参考点，远参考点的布设位置与常规远参考大地电磁法中的远参考点相同；根据勘探深度范围及测区大地背景电导率确定观测频率范围，并根据所需观测频率确定测点的同步观测时间长度、信号采样率；对于各个观测频率，根据时频转换计算单个频谱所需的时间域采样点数，确定时窗宽度，利用观测时长除以时窗宽度得到各观测频率对应的观测时窗个数；

步骤2、装置布设：

在每个测点处布设2道相互垂直的水平磁场测道、以及2道相互垂直的水平电场测道(测量通道)；各测点之间的两个观测方向均为水平x方向与水平y方向；需分析地下介质倾子矢量参数时，增加1道竖直磁场测道；

在每个远参考点处布设2道相互垂直的水平磁场测道，条件允许时可增加布设2道相互垂直的水平电场测道；不布设水平电场测道也可获得结果，但布设并记录相应数据后有利于获得更精确的结果；

每个测点处的电磁场测量装置与传统大地电磁法相同。

步骤3、数据采集：

本发明的数据采集采用时空阵列的方式；同步采集各个测点及远参考点处的电场及磁场分量数据，进行时频转换后，得到各观测频率对应的频域观测数据；各个观测频率的数据相互独立，处理方式相同；对其中任一观测频率，设其对应的观测时窗个数为I，对应的观测数据包括测点观测数据、远参考点观测数据、测点间水平磁场差分数据三个部分的数据；

第一部分、测点观测数据：包括测区内J个测点所有测道所记录的数据；设第j(j＝1,2,...,J)个测点的测道数为κ_j，包括2道相互垂直的水平磁场测道以及κ_j-2道水平电场和竖直磁场测道，所有测点的总测道数为根据测点观测数据构建测点信号时空数据矩阵X：

$X=\left(\begin{array}{cccc}{X}_{11}& X_{12}& \mathrm{..}& X_{1I}\\ X_{21}& X_{22}& \mathrm{..}& X_{2I}\\ .& .& .& .\\ .& .& .& .\\ X_{K1}& X_{K2}& \mathrm{..}& X_{KI}\end{array}\right)---\left(1\right)$

其中，X_ki为第k(k＝1,2,...,K)个测道第i(i＝1,2,...,I)个时窗的观测数据；

第二部分、远参考点观测数据：包括J'个远参考点所有K'个测道所记录的数据；根据远参考点观测数据构建天然电磁场信号时空数据矩阵X^r：

$X^r=\left(\begin{array}{cccc}{X}_{11}^r& X_{12}^r& \mathrm{..}& X_{1I}^r\\ X_{21}^r& X_{22}^r& \mathrm{..}& X_{2I}^r\\ .& .& .& .\\ .& .& .& .\\ X_{K^{\′}1}^r& X_{K^{\′}2}^r& \mathrm{..}& X_{K^{\′}I}^r\end{array}\right)---\left(2\right)$

其中，为第k(k＝1,2,...,K')个远参考测道第i(i＝1,2,...,I)个时窗的观测数据；

第三部分、测点间水平磁场差分数据：

以两测点间同一方向的水平磁场差值作为差分磁场测道；因每个测点均含有相互垂直的两个方向的水平磁场，则对J个测点，在每个方向上进行空间两两差分，共可获得差分磁场测道K＝J(J-1)个；

根据各差分磁场测道所记录的数据，即测点间水平磁场差分数据构建相关噪声场信号时空数据矩阵Y：

$Y=\left(\begin{array}{cccc}{Y}_{11}& Y_{12}& \mathrm{..}& Y_{1I}\\ Y_{21}& Y_{22}& \mathrm{..}& Y_{2I}\\ .& .& .& .\\ .& .& .& .\\ Y_{\underset{\‾}{K}1}& Y_{\underset{\‾}{K}2}& \mathrm{..}& Y_{\underset{\‾}{K}I}\end{array}\right).---\left(3\right)$

其中，Y_ki为第k(k＝1,2,...,K)个差分磁场测道第i(i＝1,2,...,I)个时窗的观测数据；

步骤4、数据处理，包括以下步骤：

第一步，对天然电磁场信号时空数据矩阵X^r进行矩阵分解求解天然电磁场源极化参数α：

设天然场源的个数为M，将X^r写为空间矩阵与时间矩阵的乘积形式：

X^r＝U^rα+ε^r,(4)

其中各参数展开式为：

$\mathrm{\α}=\left(\begin{array}{cccc}{\mathrm{\α}}_{11}& {\mathrm{\α}}_{12}& \mathrm{..}& {\mathrm{\α}}_{1I}\\ {\mathrm{\α}}_{21}& {\mathrm{\α}}_{22}& \mathrm{..}& {\mathrm{\α}}_{2I}\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\α}}_{M1}& {\mathrm{\α}}_{M2}& \mathrm{..}& {\mathrm{\α}}_{MI}\end{array}\right),U^r=\left(\begin{array}{cccc}{U}_{11}^r& U_{12}^r& \mathrm{..}& U_{1M}^r\\ U_{21}^r& U_{22}^r& \mathrm{..}& U_{2M}^r\\ .& .& .& .\\ .& .& .& .\\ U_{K^{\′}1}^r& U_{K^{\′}2}^r& \mathrm{..}& U_{K^{\′}M}^r\end{array}\right)---\left(5\right)$

其中，M×I阶矩阵α为天然场源极化参数，仅与天然场源个数及观测时窗有关，α_mi(m＝1,2,...,M；i＝1,2,...,I)为其元素，表示第m个天然场源第i个时窗的极化参数，；K'×M阶矩阵U^r为X^r对应于天然场源极化参数α的空间模数，仅与天然场源个数及空间各测道有关，为其元素，表示第k个测道对应于第m个天然场源的空间模数；ε^r为X^r中的不相关噪声；

通过对X^r进行矩阵分解获得α的估计值；

第二步，通过对相关噪声场信号时空数据矩阵Y进行矩阵分解求解相关噪声场源极化参数β：

设相关噪声场源的个数为N，将Y写为空间矩阵与时间矩阵的乘积形式：

Y＝V^cβ+ε^Y,(6)

其中各参数展开式为：

$\mathrm{\β}=\left(\begin{array}{cccc}{\mathrm{\β}}_{11}& {\mathrm{\β}}_{12}& \mathrm{..}& {\mathrm{\β}}_{1I}\\ {\mathrm{\β}}_{21}& {\mathrm{\β}}_{22}& \mathrm{..}& {\mathrm{\β}}_{2I}\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\β}}_{N1}& {\mathrm{\β}}_{N2}& \mathrm{..}& {\mathrm{\β}}_{NI}\end{array}\right),V^c=\left(\begin{array}{cccc}{V}_{11}^c& V_{12}^c& \mathrm{..}& V_{1N}^c\\ V_{21}^c& V_{22}^c& \mathrm{..}& V_{2N}^c\\ .& .& .& .\\ .& .& .& .\\ V_{\underset{\‾}{K}1}^c& V_{\underset{\‾}{K}2}^c& \mathrm{..}& V_{\underset{\‾}{K}N}^c\end{array}\right)---\left(7\right)$

其中，N×I阶矩阵β为相关噪声场源极化参数，仅与相关噪声场源个数及观测时窗有关，β_ni(n＝1,2,...,N；i＝1,2,...,I)为其元素，表示第n个相关噪声场源第i个时窗的极化参数；K×N阶矩阵V^c为Y对应于相关噪声场源极化参数β的空间模数，仅与相关噪声场源个数及空间各测道有关，为其元素，表示第k个测道对应于第n个相关噪声场源的空间模数；ε^Y为Y中的不相关噪声；

通过对Y进行矩阵分解获得β的估计值；

第三步，利用已获取的α、β估计值以及测点信号时空数据矩阵X，求解测点对应于天然场源的空间模数U和对应于相关噪声场源的空间模数V：

测点观测数据由天然场响应和相关噪声场响应两部分叠加组成，将X写为空间矩阵与时间矩阵的乘积形式，

$X=U\mathrm{\α}+V\mathrm{\β}+\mathrm{\ϵ}=\left[\begin{array}{cc}U& V\end{array}\right]\left[\begin{array}{c}\mathrm{\α}\\ \mathrm{\β}\end{array}\right]+\mathrm{\ϵ}=W\mathrm{\γ}+\mathrm{\ϵ},---\left(8\right)$

其中，W＝[UV]为总空间模数，γ＝[αβ]^*为总场源极化参数，ε为X中的不相关噪声；γ表示为：

$\mathrm{\γ}={\left[\begin{array}{cc}\mathrm{\α}& \mathrm{\β}\end{array}\right]}^{*}=\left(\begin{array}{cccc}{\mathrm{\α}}_{11}& {\mathrm{\α}}_{11}& \mathrm{..}& {\mathrm{\α}}_{1I}\\ {\mathrm{\α}}_{21}& {\mathrm{\α}}_{22}& \mathrm{..}& {\mathrm{\α}}_{2I}\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\α}}_{M1}& {\mathrm{\α}}_{M2}& \mathrm{..}& {\mathrm{\α}}_{MI}\\ {\mathrm{\β}}_{11}& {\mathrm{\β}}_{12}& \mathrm{..}& {\mathrm{\β}}_{1I}\\ {\mathrm{\β}}_{21}& {\mathrm{\β}}_{22}& \mathrm{..}& {\mathrm{\β}}_{2I}\\ .& .& .& .\\ .& .& .& .\\ {\mathrm{\β}}_{N1}& {\mathrm{\β}}_{N2}& \mathrm{..}& {\mathrm{\β}}_{NI}\end{array}\right)---\left(9\right)$

由(8)式求得W以及U、V的估计值为：

[UV]＝W(11)

上角标表示共轭转置矩阵，上角标^-1表示矩阵的逆；W的展开式为：

$W=\[\begin{array}{cc}U& V\end{array}\]=\left(\begin{array}{cccccccc}{U}_{11}& U_{12}& \mathrm{..}& U_{1M}& V_{11}& V_{12}& \mathrm{..}& V_{1N}\\ U_{21}& U_{22}& \mathrm{..}& U_{2M}& V_{21}& V_{22}& \mathrm{..}& V_{2N}\\ .& .& & .& .& .& & .\\ .& .& \mathrm{..}& .& .& .& \mathrm{..}& .\\ U_{K1}& U_{K2}& \mathrm{..}& U_{KM}& V_{K1}& V_{K2}& \mathrm{..}& V_{KN}\end{array}\right).---\left(12\right)$

$U=\left(\begin{array}{cccc}{U}_{11}& U_{12}& \mathrm{..}& U_{1M}\\ U_{21}& U_{22}& \mathrm{..}& U_{2M}\\ .& .& & .\\ .& .& \mathrm{..}& .\\ U_{K1}& U_{K2}& \mathrm{..}& U_{KM}\end{array}\right),V=\left(\begin{array}{cccc}{V}_{11}& V_{12}& \mathrm{..}& V_{1N}\\ V_{21}& V_{22}& \mathrm{..}& V_{2N}\\ .& .& & .\\ .& .& \mathrm{..}& .\\ V_{K1}& V_{K2}& \mathrm{..}& V_{KN}\end{array}\right)$

其中，U_km(k＝1,2,...,K；m＝1,2,...,M)为U的元素，表示第k个测道对应于第m个天然场源的空间模数；V_kn(k＝1,2,...,K；n＝1,2,...,N)为V的元素；

对第j(j＝1,2,...,J)个测点，设其各测道对应的标号分别为k+1,k+2,...,k+κ_j(k＝κ₁+...+κ_j-1)，则在空间模数矩阵U和V中，与第j(j＝1,2,...,J)个测点对应的空间模数U_(j)、V_(j)为：

$U_{\left(j\right)}=\left(\begin{array}{cccc}{U}_{(k+1)1}& U_{(k+1)2}& \mathrm{..}& U_{(k+1)M}\\ U_{(k+2)1}& U_{(k+2)2}& \mathrm{..}& U_{(k+2)M}\\ \begin{array}{c}.\\ .\end{array}& \begin{array}{c}.\\ .\end{array}& \mathrm{..}& \begin{array}{c}.\\ .\end{array}\\ U_{(k+{\mathrm{\κ}}_j)1}& U_{(k+{\mathrm{\κ}}_j)2}& \mathrm{..}& U_{(k+{\mathrm{\κ}}_j)M}\end{array}\right),V_{\left(j\right)}=\left(\begin{array}{cccc}{V}_{(k+1)1}& V_{(k+1)2}& \mathrm{..}& V_{(k+1)N}\\ V_{(k+2)1}& V_{(k+2)2}& \mathrm{..}& V_{(k+2)N}\\ \begin{array}{c}.\\ .\end{array}& \begin{array}{c}.\\ .\end{array}& \mathrm{..}& \begin{array}{c}.\\ .\end{array}\\ V_{(k+{\mathrm{\κ}}_j)1}& V_{(k+{\mathrm{\κ}}_j)2}& \mathrm{..}& V_{(k+{\mathrm{\κ}}_j)N}\end{array}\right)---\left(13\right)$

最后，利用求得的空间模数U和V计算各个测点的天然场阻抗张量与相关噪声场阻抗张量；对第j(j＝1,2,...,J)个测点，记天然场阻抗张量为相关噪声场阻抗张量为则

其中，分别为空间模数U_(j)、V_(j)中对应于水平磁场测道(输入端)的部分，分别为空间模数U_(j)、V_(j)中对应于水平电场测道及竖直磁场测道(输出端)的部分；具体的，对于(1)式所示的测点信号时空数据矩阵，

$U_{\left(j\right)}^H=\left(\begin{array}{cccc}{U}_{(k+1)1}& U_{(k+1)2}& \mathrm{..}& U_{(k+1)M}\\ U_{(k+2)1}& U_{(k+2)2}& \mathrm{..}& U_{(k+2)M}\end{array}\right),V_{\left(j\right)}^H=\left(\begin{array}{cccc}{V}_{(k+1)1}& V_{(k+1)2}& \mathrm{..}& V_{(k+1)N}\\ V_{(k+2)1}& V_{(k+2)2}& \mathrm{..}& V_{(k+2)N}\end{array}\right),---\left(16\right)$

$U_{\left(j\right)}^E=\left(\begin{array}{cccc}{U}_{(k+3)1}& U_{(k+3)2}& \mathrm{..}& U_{(k+3)M}\\ U_{(k+4)1}& U_{(k+4)2}& \mathrm{..}& U_{(k+4)M}\\ .& .& .& .\\ .& .& .& .\\ U_{(k+{\mathrm{\κ}}_j)1}& U_{(k+{\mathrm{\κ}}_j)2}& \mathrm{..}& U_{(k+{\mathrm{\κ}}_j)M}\end{array}\right),V_{\left(j\right)}^E=\left(\begin{array}{cccc}{V}_{(k+3)1}& V_{(k+3)2}& \mathrm{..}& V_{(k+3)N}\\ V_{(k+4)1}& V_{(k+4)2}& \mathrm{..}& V_{(k+4)N}\\ .& .& .& .\\ .& .& .& .\\ V_{(k+{\mathrm{\κ}}_j)1}& V_{(k+{\mathrm{\κ}}_j)2}& \mathrm{..}& V_{(k+{\mathrm{\κ}}_j)N}\end{array}\right).---\left(17\right)$

需说明，此处为与前述某单一观测频率数据所得的计算结果，对应了地下某一深度处的电性信息；对于不同的观测频率，可以分别独立得到相对应的进而通过接收不同频率的电磁信号，可以获得地下不同深度的电性介质分布。

与常规大地电磁法所得的天然场阻抗类似，可用于分析地下介质电性信息；但已不包含相关噪声的影响。为本方法所求得的相关噪声场信息，可用于分析测区内的信噪环境。

进一步地，所述步骤3的第一部分中，测点信号时空数据矩阵X的构建方法为：

$X=\left(\begin{array}{cccc}{X}_{11}& X_{12}& \mathrm{..}& X_{1I}\\ X_{21}& X_{22}& \mathrm{..}& X_{2I}\\ .& .& .& .\\ .& .& .& .\\ X_{K1}& X_{K2}& \mathrm{..}& X_{KI}\end{array}\right)=\left(\begin{array}{cccc}{h}_{11}& h_{12}& \mathrm{..}& h_{1I}\\ e_{11}& e_{12}& \mathrm{..}& e_{1I}\\ .& .& .& .\\ .& .& .& .\\ h_{J1}& h_{J2}& \mathrm{..}& h_{JI}\\ e_{J1}& e_{J2}& \mathrm{..}& e_{JI}\end{array}\right)---\left(18\right)$

其中，

$h_{ji}=\left(\begin{array}{c}{\left(H_{ji}\right)}_x\\ {\left(H_{ji}\right)}_y\end{array}\right),$

h_ji为第j(j＝1,2,...,J)个测点第i(i＝1,2,...,I)个时窗的水平磁场观测矢量，(H_ji)_x为其水平x方向分量，(H_ji)_y为其水平y方向分量；

当观测竖直磁场测道时，

$e_{ji}=\left(\begin{array}{c}{\left(E_{ji}\right)}_x\\ {\left(E_{ji}\right)}_y\\ {\left(H_{ji}\right)}_z\end{array}\right)$

e_ji为第j个测点第i个时窗的输出端观测矢量，(E_ji)_x为其水平x方向电场分量，(E_ji)_y为其水平y方向分量，(H_ji)_z为其竖直磁场分量；

当不观测竖直磁场测道时，

$e_{ji}=\left(\begin{array}{c}{\left(E_{ji}\right)}_x\\ {\left(E_{ji}\right)}_y\end{array}\right).$

进一步地，所述步骤3的第二部分中，天然电磁场信号时空数据矩阵X^r的构建方法为：

当远参考点处仅布设并记录2道相互垂直的水平磁场测道时，

$X^r=\left(\begin{array}{cccc}{X}_{11}^r& X_{12}^r& \mathrm{..}& X_{1I}^r\\ X_{21}^r& X_{22}^r& \mathrm{..}& X_{2I}^r\\ .& .& .& .\\ .& .& .& .\\ X_{K^{\′}1}^r& X_{K^{\′}2}^r& \mathrm{..}& X_{K^{\′}I}^r\end{array}\right)=\left(\begin{array}{cccc}{h}_{11}^r& h_{12}^r& \mathrm{..}& h_{1I}^r\\ h_{21}^r& h_{22}^r& \mathrm{..}& h_{2I}^r\\ .& .& .& .\\ .& .& .& .\\ h_{J^{\′}1}^r& h_{J^{\′}2}^r& \mathrm{..}& h_{J^{\′}I}^r\end{array}\right)---\left(19\right)$

$h_{ji}^r=\left(\begin{array}{c}{\left(H_{ji}^r\right)}_x\\ {\left(H_{ji}^r\right)}_y\end{array}\right),$

为第j(j＝1,2,...,J')个远参考点第i(i＝1,2,...,I)个时窗的水平磁场观测矢量，为其水平x方向分量，为其水平y方向分量；

当远参考点处增加布设并记录2道相互垂直的水平电场测道时，

$X^r=\left(\begin{array}{cccc}{X}_{11}^r& X_{12}^r& \mathrm{..}& X_{1I}^r\\ X_{21}^r& X_{22}^r& \mathrm{..}& X_{2I}^r\\ .& .& .& .\\ .& .& .& .\\ X_{K^{\′}1}^r& X_{K^{\′}2}^r& \mathrm{..}& X_{K^{\′}I}^r\end{array}\right)=\left(\begin{array}{cccc}{h}_{11}^r& h_{12}^r& \mathrm{..}& h_{1I}^r\\ e_{11}^r& e_{12}^r& \mathrm{..}& e_{1I}^r\\ .& .& .& .\\ .& .& .& .\\ h_{J^{\′}1}^r& h_{J^{\′}2}^r& \mathrm{..}& h_{J^{\′}I}^r\\ e_{J^{\′}1}^r& e_{J^{\′}2}^r& \mathrm{..}& e_{J^{\′}I}^r\end{array}\right)---\left(20\right)$

其中

$e_{ji}^r=\left(\begin{array}{c}{\left(E_{ji}^r\right)}_x\\ {\left(E_{ji}^r\right)}_y\end{array}\right)$

为第j(j＝1,2,...,J')个远参考点第i(i＝1,2,...,I)个时窗的水平电场观测矢量，为其水平x方向分量，为其水平y方向分量。

进一步地，所述步骤3的第三部分中，设第k(k＝1,2,...,K)个差分磁场测道由第j和j'(j,j'＝1,2,...,J且j≠j')两个测点的d方向(d方向为x、y两个方向之一)水平磁场进行差分所得，则第k个差分磁场测道第i个时窗的水平磁场差分数据Y_ki为：

Y_ki＝H_j,d,i-H_j',d,i,(21)

其中，H_j,d,i、H_j',d,i分别为第j和j'个测点的d方向第i个时窗的水平磁场观测数据。

进一步地，所述步骤4的第一步和第二步中，矩阵分解方法可采用奇异值分解、主成分分析法和稳健性主成分分析法等。为便于说明，此处以形式最简单的奇异值分解为例；具体步骤为：

[U,S,V]＝svd(X^r),(22)

U^r＝(U)^K'×M,α＝(SV ^*)^M×I,(23)

其中，svd表示矩阵的奇异值分解，U,S,V分别为奇异值分解所得的参数矩阵，上角标^*表示矩阵转置；(U)^K'×M表示U的1～M列组成的K'×M阶矩阵，(SV ^*)^M×I表示SV ^*的1～M行组成的M×I阶矩阵。

进一步地，所述步骤4的第二步中，采用的矩阵分解方法为奇异值分解；具体步骤为：

[U,S,V]＝svd(Y),(24)

V^c＝(U) ^K×N,β＝(SV ^*)^N×I,(25)

其中，svd表示矩阵的奇异值分解，U,S,V分别为奇异值分解所得的参数矩阵，上角标^*表示矩阵转置；(U) ^K×N表示U的1～N列组成的K×N阶矩阵，(SV ^*)^N×I表示SV ^*的1～N行组成的N×I阶矩阵。

本发明所述的时空阵列电磁方法理论的具体过程为：

对阵列观测数据，假设阵列数据中共含测点J个，且J≥2，第j(j＝1,2,...,J)个测点的测道数为κ_j，包括2道相互垂直的水平磁场测道以及κ_j-2道水平电场和竖直磁场测道，总测道数为共有同步观测时窗I个，所有测道、时窗中均允许存在噪声。

在考虑相关噪声的条件下，假设观测系统中共有L个场源，包含M个天然电磁场源以及N个相关噪声场源，L＝M+N；在传统大地电磁法中，通常假设天然场为均匀平面波电磁场，此时可简化认为M＝2。在本发明中，同样认为M＝2。天然电磁场源随时间变化的极化因子为α，相关噪声场源随时间变化的极化因子为β。

电磁勘探的模型为线性时不变系统，涉及的场源、大地以及观测点等对象均满足线性时不变系统的基本特征。

具体的，首先，时不变系统的系统参数不随时间而变化，或者说激励与响应的时间变化一致。因此，对某一天然电磁场源α_m(m＝1,2,...,M)，观测响应满足

$X_{ki,m}^{MT}=U_{km}{\mathrm{\α}}_{mi},---\left(26\right)$

其中，i＝1,2,...,I为观测时窗序号，m＝1,2,...,M为天然场源标号，k＝1,2,...,K为电场、磁场测道编号，α_mi为第m个天然场源第i个时窗的极化参数，表征了场源的时变信息，不随测点的空间位置变化；U_km为第k个测道对应于第m个天然场源的空间模数，是大地电性参数ρ、测点与第m个场源的距离r以及观测角度等参数的函数，不随时间变化；为第k个测道第i个时窗对应于第m个天然场源的天然场响应。

同样的，对某一相关噪声场源β_n(n＝1,2,...,N)，观测响应满足

$X_{ki,n}^{CN}=V_{kn}{\mathrm{\β}}_{ni},---\left(27\right)$

其中，n＝1,2,...,N为相关噪声场源的标号，β_ni为第n个相关噪声场源第i个时窗的极化参数，不随测点的空间位置变化；V_kn为第k个测道对应于第n个相关噪声场源的空间模数，不随时间变化；为第k个测道第i个时窗对应于第n个相关噪声场源的噪声响应。

其次，线性系统满足叠加原理，不同场源的综合作用在观测点处可视为各个单一场源作用的叠加。因此，实际观测响应满足

$X_{ki}=\underset{m=1}{\overset{M}{\Σ}}{X}_{ki,m}^{MT}+\underset{n=1}{\overset{N}{\Σ}}{X}_{ki,n}^{CN}+{\mathrm{\ϵ}}_{ki}=\underset{m=1}{\overset{M}{\Σ}}{U}_{km}{\mathrm{\α}}_{mi}+\underset{n=1}{\overset{N}{\Σ}}{V}_{kn}{\mathrm{\β}}_{ni}+{\mathrm{\ϵ}}_{ki},---\left(28\right)$

其中，ε_ki为不相关噪声，X_ki为第k(k＝1,2,...,K)个测道第i(i＝1,2,...,I)个时窗的实际观测响应。

(28)式可写为矩阵形式，

X＝Uα+Vβ+ε,(29)

此即技术方案中的(8)式；其中，X为测点信号时空数据矩阵，由测点的观测磁场和电场组成，α为天然场源极化参数，β为相关噪声场源极化参数；U为对应于天然场源的空间模数，V为对应于相关噪声场源的空间模数，ε为不相关噪声矩阵。各参数的展开式如(1)、(5)、(7)式所示。

在远参考点处，测区内的相关噪声场已没有影响，观测响应由天然场响应和不相关噪声组成。与(29)式的推导类似，可以将天然电磁场信号时空数据矩阵X^r写成时空矩阵相乘的形式，

X^r＝U^rα+ε^r(30)

此即技术方案中的(4)式；其中，U^r为X^r对应于天然场源极化参数α的空间模数，仅与天然场源个数及空间各测道有关；ε^r为X^r中的不相关噪声。参数的展开式如(2)、(5)式所示。

在测区内，含噪测点j的磁场满足：

其中，表示测点j的观测磁场矢量(包含两个相互垂直的水平分量)，表示其中的天然电磁场成分，表示相关噪声成分，表示不相关噪声成分。

由于天然电磁场源距离测区足够远，其输入信号(水平磁场)在相当空间范围内可视为均匀平面电磁波，因此我们可对水平天然磁场作如下假设：

其中j、j'分别表示两个不同的测点，j≠j'。由(31)、(32)式两测点j、j'水平磁场的差分信号可写为，

其中，表示的相关噪声部分，表示的不相关噪声部分。

不难发现(33)式中仅含有相关噪声与不相关噪声；进一步的，将所有含噪测点的水平磁场两两相减可构成仅包含相关噪声项的时空阵列数据集Y；其形式如(3)式所示。与(29)、(30)两式相同，相关噪声场信号时空数据矩阵Y可以写成时空矩阵相乘的形式，

Y＝V^cβ+ε^Y,(34)

此即技术方案中的(6)式；其中，V^c为Y对应于相关噪声场极化参数β的空间模数，ε^Y为不相关噪声项。参数的展开式如(7)式所示。

(29)、(30)及(34)式即为本方法所述的时空阵列电磁信噪统计模型。

频域率电磁勘探方法关心的是每个观测点处不随时间变化的天然场空间模数U，并从U中提取地下电性参数ρ。从时空阵列电磁信噪统计模型中提取地电参数的求解策略可分为四个步骤。

第一步，利用天然电磁场信号时空数据矩阵X^r获得天然场源的时间变化项α。对于形如(30)式所示的时空乘积矩阵，可以采用矩阵分解的方法进行求解，一种合理的矩阵分解方法是稳健的主成分分析方法。在本文中，矩阵分解方法以奇异值分解为例，(30)式的求解方法为

[U,S,V]＝svd(X^r),(35)

其中，svd表示矩阵的奇异值分解，U,S,V分别为奇异值分解所得的参数矩阵，上角标*表示矩阵转置。具体的，对K'×I阶矩阵Xr，奇异值分解所得参数矩阵的展开式为，

$\underset{\‾}{U}=\left(\begin{array}{cccc}{\underset{\‾}{U}}_{11}& {\underset{\‾}{U}}_{12}& \mathrm{..}& {\underset{\‾}{U}}_{1K^{\′}}\\ {\underset{\‾}{U}}_{21}& {\underset{\‾}{U}}_{22}& \mathrm{..}& {\underset{\‾}{U}}_{2K^{\′}}\\ .& .& .& .\\ .& .& .& .\\ {\underset{\‾}{U}}_{K^{\′}1}& {\underset{\‾}{U}}_{K^{\′}2}& \mathrm{..}& {\underset{\‾}{U}}_{K^{\′}{K}^{\′}}\end{array}\right),\underset{\‾}{V}=\left(\begin{array}{cccc}{\underset{\‾}{V}}_{11}& {\underset{\‾}{V}}_{12}& \mathrm{..}& {\underset{\‾}{V}}_{1I}\\ {\underset{\‾}{V}}_{21}& {\underset{\‾}{V}}_{22}& \mathrm{..}& {\underset{\‾}{V}}_{2I}\\ .& .& .& .\\ .& .& .& .\\ {\underset{\‾}{V}}_{I1}& {\underset{\‾}{V}}_{I2}& \mathrm{..}& {\underset{\‾}{V}}_{II}\end{array}\right),---\left(36\right)$

当K'<I时，K'×I阶矩阵S为

$\underset{\&OverBar;}{S}=\left(\begin{array}{cccccc}{\underset{\&OverBar;}{S}}_{11}& 0& \mathrm{..}& 0& 0& \mathrm{..}\\ 0& {\underset{\&OverBar;}{S}}_{22}& \mathrm{..}& 0& 0& \mathrm{..}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& 0& \mathrm{..}& {\underset{\&OverBar;}{S}}_{K^{\&prime;}{K}^{\&prime;}}& 0& \mathrm{..}\end{array}\right),---\left(37\right)$

当K'≥I时，K'×I阶矩阵S为

$\left(\begin{array}{cccc}{\underset{\&OverBar;}{S}}_{11}& 0& \mathrm{..}& 0\\ 0& {\underset{\&OverBar;}{S}}_{22}& \mathrm{..}& 0\\ .& .& .& .\\ .& .& .& .\\ 0& 0& \mathrm{..}& {\underset{\&OverBar;}{S}}_{II}\\ 0& 0& \mathrm{..}& 0\\ .& .& .& .\\ .& .& .& .\end{array}\right).---\left(38\right)$

于是可求得

U^r＝(U)^M,α＝(SV ^*)^M,(39)(U)^M表示U的1～M列，(SV ^*)^M表示SV ^*的1～M行；展开式为，

$U^r={\left(\underset{\&OverBar;}{U}\right)}^M=\left(\begin{array}{cccc}{\underset{\&OverBar;}{U}}_{11}& {\underset{\&OverBar;}{U}}_{12}& \mathrm{..}& {\underset{\&OverBar;}{U}}_{1M}\\ {\underset{\&OverBar;}{U}}_{21}& {\underset{\&OverBar;}{U}}_{22}& \mathrm{..}& {\underset{\&OverBar;}{U}}_{2M}\\ .& .& .& .\\ .& .& .& .\\ {\underset{\&OverBar;}{U}}_{K^{\&prime;}1}& {\underset{\&OverBar;}{U}}_{K^{\&prime;}2}& \mathrm{..}& {\underset{\&OverBar;}{U}}_{K^{\&prime;}M}\end{array}\right),---\left(40\right)$

$\mathrm{\&alpha;}={\left(\underset{\&OverBar;}{S}{\underset{\&OverBar;}{V}}^{*}\right)}^M=\left(\begin{array}{cccc}{\underset{\&OverBar;}{S}}_{11}& 0& \mathrm{..}& 0\\ 0& {\underset{\&OverBar;}{S}}_{22}& \mathrm{..}& 0\\ .& .& .& .\\ .& .& .& .\\ 0& 0& \mathrm{..}& {\underset{\&OverBar;}{S}}_{MM}\end{array}\right)\left(\begin{array}{cccc}{\underset{\&OverBar;}{V}}_{11}& {\underset{\&OverBar;}{V}}_{12}& \mathrm{..}& {\underset{\&OverBar;}{V}}_{1I}\\ {\underset{\&OverBar;}{V}}_{21}& {\underset{\&OverBar;}{V}}_{22}& \mathrm{..}& {\underset{\&OverBar;}{V}}_{2I}\\ .& .& .& .\\ .& .& .& .\\ {\underset{\&OverBar;}{V}}_{M1}& {\underset{\&OverBar;}{V}}_{M2}& \mathrm{..}& {\underset{\&OverBar;}{V}}_{MI}\end{array}\right).---\left(41\right)$

第二步，利用相关噪声场信号时空数据矩阵Y以及(34)式所示的相关噪声时空阵列统计模型，求解相关噪声场源极化参数β。与第一步类似，利用矩阵分解的方法获得对β的估计，以奇异值分解为例，

[U,S,V]＝svd(Y),(42)

其中，svd表示矩阵的奇异值分解，U,S,V分别为奇异值分解所得的参数矩阵，上角标*表示矩阵转置。具体的，对K×I阶矩阵Y，奇异值分解所得参数矩阵的展开式为，

$\underset{\&OverBar;}{U}=\left(\begin{array}{cccc}{\underset{\&OverBar;}{U}}_{11}& {\underset{\&OverBar;}{U}}_{12}& \mathrm{..}& {\underset{\&OverBar;}{U}}_{1\underset{\&OverBar;}{K}}\\ {\underset{\&OverBar;}{U}}_{21}& {\underset{\&OverBar;}{U}}_{22}& \mathrm{..}& {\underset{\&OverBar;}{U}}_{2\underset{\&OverBar;}{K}}\\ .& .& .& .\\ .& .& .& .\\ {\underset{\&OverBar;}{U}}_{\underset{\&OverBar;}{K}1}& {\underset{\&OverBar;}{U}}_{\underset{\&OverBar;}{K}2}& \mathrm{..}& {\underset{\&OverBar;}{U}}_{\underset{\&OverBar;}{K}\underset{\&OverBar;}{K}}\end{array}\right),\underset{\&OverBar;}{V}=\left(\begin{array}{cccc}{\underset{\&OverBar;}{V}}_{11}& {\underset{\&OverBar;}{V}}_{12}& \mathrm{..}& {\underset{\&OverBar;}{V}}_{1I}\\ {\underset{\&OverBar;}{V}}_{21}& {\underset{\&OverBar;}{V}}_{22}& \mathrm{..}& {\underset{\&OverBar;}{V}}_{2I}\\ .& .& .& .\\ .& .& .& .\\ {\underset{\&OverBar;}{V}}_{I1}& {\underset{\&OverBar;}{V}}_{I2}& \mathrm{..}& {\underset{\&OverBar;}{V}}_{II}\end{array}\right),---\left(43\right)$

当K<I时，K×I阶矩阵S为

$\underset{\&OverBar;}{S}=\left(\begin{array}{cccccc}{\underset{\&OverBar;}{S}}_{11}& 0& \mathrm{..}& 0& 0& \mathrm{..}\\ 0& {\underset{\&OverBar;}{S}}_{22}& \mathrm{..}& 0& 0& \mathrm{..}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& 0& \mathrm{..}& {\underset{\&OverBar;}{S}}_{\underset{\&OverBar;}{K}\underset{\&OverBar;}{K}}& 0& \mathrm{..}\end{array}\right),---\left(44\right)$

当K≥I时，K×I阶矩阵S为

$\left(\begin{array}{cccc}{\underset{\&OverBar;}{S}}_{11}& 0& \mathrm{..}& 0\\ 0& {\underset{\&OverBar;}{S}}_{22}& \mathrm{..}& 0\\ .& .& .& .\\ .& .& .& .\\ 0& 0& \mathrm{..}& {\underset{\&OverBar;}{S}}_{II}\\ 0& 0& \mathrm{..}& 0\\ .& .& .& .\\ .& .& .& .\end{array}\right).---\left(45\right)$

于是可求得

V^c＝(U)^N,β＝(SV ^*)^N,(46)

(U)^N表示U的1～N列，(SV ^*)^N表示SV ^*的1～N行；展开式为，

$V^c={\left(\underset{\&OverBar;}{U}\right)}^N=\left(\begin{array}{cccc}{\underset{\&OverBar;}{U}}_{11}& {\underset{\&OverBar;}{U}}_{12}& \mathrm{..}& {\underset{\&OverBar;}{U}}_{1N}\\ {\underset{\&OverBar;}{U}}_{21}& {\underset{\&OverBar;}{U}}_{22}& \mathrm{..}& {\underset{\&OverBar;}{U}}_{2N}\\ .& .& .& .\\ .& .& .& .\\ {\underset{\&OverBar;}{U}}_{\underset{\&OverBar;}{K}1}& {\underset{\&OverBar;}{U}}_{\underset{\&OverBar;}{K}2}& \mathrm{..}& {\underset{\&OverBar;}{U}}_{\underset{\&OverBar;}{K}N}\end{array}\right),---\left(47\right)$

$\mathrm{\&beta;}={\left(\underset{\&OverBar;}{S}{\underset{\&OverBar;}{V}}^{*}\right)}^N=\left(\begin{array}{cccc}{\underset{\&OverBar;}{S}}_{11}& 0& \mathrm{..}& 0\\ 0& {\underset{\&OverBar;}{S}}_{22}& \mathrm{..}& 0\\ .& .& .& .\\ .& .& .& .\\ 0& 0& \mathrm{..}& {\underset{\&OverBar;}{S}}_{NN}\end{array}\right)\left(\begin{array}{cccc}{\underset{\&OverBar;}{V}}_{11}& {\underset{\&OverBar;}{V}}_{12}& \mathrm{..}& {\underset{\&OverBar;}{V}}_{1I}\\ {\underset{\&OverBar;}{V}}_{21}& {\underset{\&OverBar;}{V}}_{22}& \mathrm{..}& {\underset{\&OverBar;}{V}}_{2I}\\ .& .& .& .\\ .& .& .& .\\ {\underset{\&OverBar;}{V}}_{N1}& {\underset{\&OverBar;}{V}}_{N2}& \mathrm{..}& {\underset{\&OverBar;}{V}}_{NI}\end{array}\right).---\left(48\right)$

第三步，利用前两步所求得的天然场源极化参数α以及相关噪声场源极化参数β，通过测点信噪统计模型(29)式估计U和V。(29)式可写为，

$X=\left[\begin{array}{cc}U& V\end{array}\right]\left[\begin{array}{c}\mathrm{\&alpha;}\\ \mathrm{\&beta;}\end{array}\right]+\mathrm{\&epsiv;}=W\mathrm{\&gamma;}+\mathrm{\&epsiv;},---\left(49\right)$

W的估计值为，

[UV]＝W,(51)

上角标表示共轭转置矩阵，上角标^-1表示矩阵的逆。各参数的展开式如(1)、(9)、(12)式所示。

最后，利用求得的空间模数U和V计算天然场平面波阻抗与相关噪声场非平面波阻抗。对第j(j＝1,2,...,J)个测点，假设其对应的测道标号分别为k+1,k+2,...,k+κ_j,(k＝κ₁+...+κ_j-1)。于是在空间模数矩阵U和V中，与第j(j＝1,2,...,J)个测点对应的空间模数U_(j)、V_(j)为

$U_{\left(j\right)}=\left(\begin{array}{cccc}{U}_{(k+1)1}& U_{(k+1)2}& \mathrm{..}& U_{(k+1)M}\\ U_{(k+2)1}& U_{(k+2)2}& \mathrm{..}& U_{(k+2)M}\\ .& .& .& .\\ .& .& .& .\\ U_{(k+{\mathrm{\&kappa;}}_j)1}& U_{(k+{\mathrm{\&kappa;}}_j)2}& \mathrm{..}& U_{(k+{\mathrm{\&kappa;}}_j)M}\end{array}\right),V_{\left(j\right)}=\left(\begin{array}{cccc}{V}_{(k+1)1}& V_{(k+1)2}& \mathrm{..}& V_{(k+1)N}\\ V_{(k+2)1}& V_{(k+2)2}& \mathrm{..}& V_{(k+2)N}\\ .& .& .& .\\ .& .& .& .\\ V_{(k+{\mathrm{\&kappa;}}_j)1}& V_{(k+{\mathrm{\&kappa;}}_j)2}& \mathrm{..}& V_{(k+{\mathrm{\&kappa;}}_j)N}\end{array}\right).---\left(52\right)$

空间模数U_(j)、V_(j)中对应于水平磁场测道(输入端)的部分标示为对应于水平电场测道及竖直磁场测道(输出端)的部分为标示具体的，对于(1)式所示的测点信号时空数据矩阵，

$U_{\left(j\right)}^H=\left(\begin{array}{cccc}{U}_{(k+1)1}& U_{(k+1)2}& \mathrm{..}& U_{(k+1)M}\\ U_{(k+2)1}& U_{(k+2)2}& \mathrm{..}& U_{(k+2)M}\end{array}\right),V_{\left(j\right)}^H=\left(\begin{array}{cccc}{V}_{(k+1)1}& V_{(k+1)2}& \mathrm{..}& V_{(k+1)N}\\ V_{(k+2)1}& V_{(k+2)2}& \mathrm{..}& V_{(k+2)N}\end{array}\right),---\left(53\right)$

$U_{\left(j\right)}^E=\left(\begin{array}{cccc}{U}_{(k+3)1}& U_{(k+3)2}& \mathrm{..}& U_{(k+3)M}\\ U_{(k+4)1}& U_{(k+4)2}& \mathrm{..}& U_{(k+4)M}\\ .& .& .& .\\ .& .& .& .\\ U_{(k+{\mathrm{\&kappa;}}_j)1}& U_{(k+{\mathrm{\&kappa;}}_j)2}& \mathrm{..}& U_{(k+{\mathrm{\&kappa;}}_j)M}\end{array}\right),V_{\left(j\right)}^E=\left(\begin{array}{cccc}{V}_{(k+3)1}& V_{(k+3)2}& \mathrm{..}& V_{(k+3)N}\\ V_{(k+4)1}& V_{(k+4)2}& \mathrm{..}& V_{(k+4)N}\\ .& .& .& .\\ .& .& .& .\\ V_{(k+{\mathrm{\&kappa;}}_j)1}& V_{(k+{\mathrm{\&kappa;}}_j)2}& \mathrm{..}& V_{(k+{\mathrm{\&kappa;}}_j)N}\end{array}\right)---\left(54\right)$

定义天然场阻抗张量使之满足

$U_{\left(j\right)}^E=Z_j^{MT}{U}_{\left(j\right)}^H,---\left(55\right)$

或

由于奇异值分解满足条件I为单位矩阵，故(56)式可简化为

同样的，可定义相关噪声场阻抗张量

即为待求解的包含地下电性信息的电磁勘探解释参数，可进一步转化为更为直观的视电阻率、相位参数：

其中，ω为角频率，μ为磁导率；Z为张量或的元素，ρ、为张量视电阻率ρ_j、相位的元素。

有益效果：

本发明提供了一种时空阵列差分电磁勘探方法，阵列数据采集的空间布设形式灵活多样；基于统一的数据方程，尽可能的挖掘时空阵列数据中的隐含信息，且所有测道、时窗均允许存在噪声，可获得所有测点的天然场平面波阻抗与相关噪声场阻抗。利用本发明，可以在受相关噪声场干扰的区域(如我国东部地区)开展天然场源电磁法勘探。

使用本发明，采用时空阵列电磁法，在时间、空间尺度上均进行阵列数据采集，并从时空阵列数据中提取有效信息，挖掘地电参数；通过接收不同频率的电磁信号，可以获得地下不同深度的电性介质分布，可以查明地下地电特性分布、地质构造及矿产分布或者解决其它工程、水文及环境地质问题。

附图说明

图1为观测系统简化信噪模型。其中，Srx、Sry为天然电磁场源，Sc为相关噪声场源，Scx、Scy分别为Sc在平面直角坐标系中的x、y分量；测点A，B，...为分布在测区内的空间观测阵列，观测响应中包含Srx、Sry、Scx及Scy的作用；RR为远参考点，距离Scx、Scy足够远，相关噪声的影响可忽略，观测响应中主要为Srx、Sry的作用。

图2为本方法对图1所示模型观测数据的阻抗视电阻率、相位估计结果示意；其中，模型为地下电阻率为ρ₀的均匀半空间，横坐标为观测频率，纵坐标为计算数据；ρ_a/ρ₀表示根据公式(59)计算的视电阻率与真实电阻率ρ₀的比值，表示根据公式(59)计算的视相位；‘A2-MT’、‘A2-CN’分别表示包含2个测点的观测阵列，采用本文方法得到的天然场阻抗结果、相关噪声阻抗结果；该示意中取阻抗张量的xy分量结果进行展示。

图3为本方法对包含4个相关噪声场源模型的阻抗视电阻率、相位估计结果示意；其中，模型为地下电阻率为ρ₀的均匀半空间，横坐标为观测频率，纵坐标为计算数据；ρ_a/ρ₀表示根据公式(59)计算的视电阻率与真实电阻率ρ₀的比值，表示根据公式(59)计算的视相位；‘A3-MT’、‘A3-CN’分别表示包含3个测点的观测阵列，采用本文方法得到的天然场阻抗结果、相关噪声阻抗结果；该示意中取阻抗张量的xy分量结果进行展示。

具体实施方式

以下结合附图和具体实施方式对本发明作进一步的说明。

本发明所提供的时空阵列差分电磁勘探方法包括以下步骤：

(1)观测设计：确定观测目标及深度范围，根据实际勘探深度需要及测区大地背景电导率确定观测频率范围，并根据所需观测频率确定测点的同步观测时间长度及信号采样率；根据探测目标及测区范围等设计测线、测点；如图1所示，在测区内部署阵列观测点A，B，...，同时在距离测区足够远的低噪区布置远参考点RR；

(2)装置布设：利用多台电磁场测量装置，在测区内布设多个测点(大于2个)；同时在距离测区足够远的低噪区，布设至少1个远参考点；每个测点处的装置与传统大地电磁法相同，布设相互垂直的2道水平磁场测道、2道水平电场测道，必要时可增加1道垂直磁场测道；

(3)数据采集：利用GPS、北斗或其他同步方法，同步采集多个测点及远参考点处的电场及磁场分量数据；

(4)数据处理：利用观测的电磁场数据，根据本发明所提供的处理步骤及公式求解不同频率的天然场张量阻抗及相关噪声场张量阻抗；

(5)后期处理：根据数据处理所获得的天然场张量阻抗进行数据分析、反演成图及资料解释。

图2为本方法对图1所示模型观测数据的阻抗视电阻率、相位估计结果示意；其中，模型地下电阻率为ρ₀，横坐标为观测频率，纵坐标为计算数据；ρ_a/ρ₀表示根据公式(59)计算的视电阻率与真实电阻率ρ₀的比值，表示根据公式(59)计算的视相位；‘A2-MT’、‘A2-CN’分别表示包含2个测点的观测阵列，采用本文方法得到的天然场阻抗结果、相关噪声阻抗结果；图中取阻抗张量的xy分量结果进行展示。不难发现，利用本文所提供的时空阵列电磁法，获得了对应于天然场张量阻抗的视电阻率、相位以及对应于相关噪声场张量阻抗的视电阻率、相位，天然场张量阻抗具有频率测深的意义，随着频率的降低，天然场张量阻抗视电阻率ρ_a与地下真实电阻率ρ₀一致(ρ_a/ρ₀＝1)，相位稳定为45°，表明所求得的天然场张量阻抗已不受相关噪声的影响，可以获得不受畸变的的地下电性信息；同时相关噪声场张量阻抗视电阻率ρ_a与地下真实电阻率ρ₀相比在低频差异巨大(ρ_a/ρ₀>1)，相位逐渐趋于0，结合人工源电磁法的已有知识可知，所求得的相关噪声场张量阻抗准确的提取到了相关噪声场源的影响，实现了天然场阻抗与相关噪声场阻抗的分离。

当场源分布较复杂时，必须采用更大的阵列才能获得主要的非平面波极化参数及天然场张量阻抗。图3为本方法对包含4个相关噪声场源模型的阻抗视电阻率、相位估计结果示意；其中，模型地下电阻率为ρ₀，横坐标为观测频率，纵坐标为计算数据；ρ_a/ρ₀表示根据公式(59)计算的视电阻率与真实电阻率ρ₀的比值，表示根据公式(59)计算的视相位；‘A3-MT’、‘A3-CN’分别表示包含3个测点的观测阵列，采用本文方法得到的天然场阻抗结果、相关噪声阻抗结果；图中取阻抗张量的xy分量结果进行展示。与图(2)类似，不难发现，利用本文所提供的时空阵列电磁法，获得了对应于天然场张量阻抗的视电阻率、相位以及对应于相关噪声场张量阻抗的视电阻率、相位；天然场张量阻抗不受相关噪声的影响，获得了不受畸变的的地下电性信息，相关噪声场张量阻抗准确的提取到了相关噪声场源的影响，实现了天然场阻抗与相关噪声场阻抗的分离。

分析表明，采用时空阵列电磁法，在时间、空间尺度上均进行阵列数据采集，并从时空阵列数据中提取有效信息，挖掘地电参数；通过接收不同频率的电磁信号，可以获得地下不同深度的电性介质分布，达到电磁勘探的目的。

Claims (7)

1.一种时空阵列差分电磁勘探方法，其特征在于，包括以下步骤：

步骤1、观测设计：确定探测目标及勘探深度范围，在测区内设置J个同步观测测点，J≥2；同时设置J'(J'≥1)个远参考点；根据勘探深度范围及测区大地背景电导率确定观测频率范围，并根据所需观测频率确定测点的同步观测时间长度及信号采样率；对于各个观测频率，根据时频转换计算单个频谱所需的时间域采样点数，确定时窗宽度，利用观测时长除以时窗宽度得到各频率对应的观测时窗个数；

步骤2、装置布设：

在每个测点处布设2道相互垂直的水平磁场测道、以及2道相互垂直的水平电场测道(测量通道)；各测点之间的两个观测方向均为水平x方向和水平y方向；需分析地下介质倾子矢量参数时，增加1道竖直磁场测道；

在每个远参考点处布设2道相互垂直的水平磁场测道，条件允许时可增加布设2道相互垂直的水平电场测道；

步骤3、数据采集：

同步采集各个测点及远参考点处的电场及磁场分量数据，进行时频转换后，得到各观测频率对应的频域观测数据；对其中任一观测频率，设其对应的观测时窗个数为I，对应的观测数据包括测点观测数据、远参考点观测数据、测点间水平磁场差分数据三个部分的数据；

第一部分、测点观测数据：包括测区内J个测点所有测道所记录的数据；设第j(j＝1,2,...,J)个测点的测道数为κ_j，所有测点的总测道数为根据测点观测数据构建测点信号时空数据矩阵X：

$X=\left(\begin{array}{cccc}{X}_{11}& X_{12}& \&CenterDot;\&CenterDot;& X_{1I}\\ X_{21}& X_{22}& \&CenterDot;\&CenterDot;& X_{2I}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ X_{K1}& X_{K2}& \&CenterDot;\&CenterDot;& X_{KI}\end{array}\right)---\left(1\right)$

其中，X_ki为第k(k＝1,2,...,K)个测道第i(i＝1,2,...,I)个时窗的观测数据；

第二部分、远参考点观测数据：包括J'个远参考点所有K'个测道所记录的数据；设第j(j＝1,2,...,J')个远参考点的测道数为κ_j'，总测道数为根据远参考点观测数据构建天然电磁场信号时空数据矩阵X^r：

$X^r=\left(\begin{array}{cccc}{X}_{\mathrm{11}}^r& X_{\mathrm{12}}^r& \&CenterDot;\&CenterDot;& X_{\mathrm{1I}}^r\\ X_{\mathrm{21}}^r& X_{22}^r& \&CenterDot;\&CenterDot;& X_{2I}^r\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ X_{K^{\&prime;}1}^r& X_{K^{\&prime;}2}^r& \&CenterDot;\&CenterDot;& X_{K^{\&prime;}I}^r\end{array}\right)---\left(2\right)$

其中，为第k(k＝1,2,...,K')个远参考测道第i(i＝1,2,...,I)个时窗的观测数据；

第三部分、测点间水平磁场差分数据：

以两测点间同一方向的水平磁场差值作为差分磁场测道；对J个测点，在每个方向上进行空间两两差分，共可获得差分磁场测道K＝J(J-1)个；

根据各差分磁场测道所记录的数据，即测点间水平磁场差分数据构建相关噪声场信号时空数据矩阵Y：

$Y=\left(\begin{array}{cccc}{Y}_{11}& Y_{12}& \&CenterDot;\&CenterDot;& Y_{1I}\\ Y_{21}& Y_{22}& \&CenterDot;\&CenterDot;& Y_{2I}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ Y_{\underset{\&OverBar;}{K}1}& Y_{\underset{\&OverBar;}{K}2}& \&CenterDot;\&CenterDot;& Y_{\underset{\&OverBar;}{K}I}\end{array}\right).---\left(3\right)$

其中，Y_ki为第k(k＝1,2,...,K)个差分磁场测道第i(i＝1,2,...,I)个时窗的观测数据；

步骤4、数据处理，包括以下步骤：

第一步，对天然电磁场信号时空数据矩阵X^r进行矩阵分解求解天然电磁场源极化参数α：

设天然场源的个数为M，将X^r写为空间矩阵与时间矩阵的乘积形式：

X^r＝U^rα+ε^r,(4)

其中各参数展开式为：

$\mathrm{\&alpha;}=\left(\begin{array}{cccc}{\mathrm{\&alpha;}}_{11}& {\mathrm{\&alpha;}}_{12}& \&CenterDot;\&CenterDot;& {\mathrm{\&alpha;}}_{1I}\\ {\mathrm{\&alpha;}}_{21}& {\mathrm{\&alpha;}}_{22}& \&CenterDot;\&CenterDot;& {\mathrm{\&alpha;}}_{2I}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ {\mathrm{\&alpha;}}_{M1}& {\mathrm{\&alpha;}}_{M2}& \&CenterDot;\&CenterDot;& {\mathrm{\&alpha;}}_{MI}\end{array}\right),U^r=\left(\begin{array}{cccc}{U}_{11}^r& U_{12}^r& \&CenterDot;\&CenterDot;& U_{1M}^r\\ U_{21}^r& U_{22}^r& \&CenterDot;\&CenterDot;& U_{2M}^r\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ U_{K^{\&prime;}1}^r& U_{K^{\&prime;}2}^r& \&CenterDot;\&CenterDot;& U_{K^{\&prime;}M}^r\end{array}\right)---\left(5\right)$

其中，M×I阶矩阵α为天然场源极化参数，仅与天然场源个数及观测时窗有关，α_mi(m＝1,2,...,M；i＝1,2,...,I)为其元素，表示第m个天然场源第i个时窗的极化参数；K'×M阶矩阵U^r为X^r对应于天然场源极化参数α的空间模数，仅与天然场源个数及空间各测道有关，为其元素，表示第k个测道对应于第m个天然场源的空间模数；ε^r为X^r中的不相关噪声；

通过对X^r进行矩阵分解获得α的估计值；

第二步，通过对相关噪声场信号时空数据矩阵Y进行矩阵分解求解相关噪声场源极化参数β：

设相关噪声场源的个数为N，将Y写为空间矩阵与时间矩阵的乘积形式：

Y＝V^cβ+ε^Y,(6)

其中各参数展开式为：

$\mathrm{\&beta;}=\left(\begin{array}{cccc}{\mathrm{\&beta;}}_{11}& {\mathrm{\&beta;}}_{12}& \&CenterDot;\&CenterDot;& {\mathrm{\&beta;}}_{1I}\\ {\mathrm{\&beta;}}_{21}& {\mathrm{\&beta;}}_{22}& \&CenterDot;\&CenterDot;& {\mathrm{\&beta;}}_{2I}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ {\mathrm{\&beta;}}_{N1}& {\mathrm{\&beta;}}_{N2}& \&CenterDot;\&CenterDot;& {\mathrm{\&beta;}}_{NI}\end{array}\right),V^c=\left(\begin{array}{cccc}{V}_{11}^c& V_{12}^c& \&CenterDot;\&CenterDot;& V_{1N}^c\\ V_{21}^c& V_{\mathrm{22}}^c& \&CenterDot;\&CenterDot;& V_{2N}^c\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ V_{\underset{\&OverBar;}{K}1}^c& V_{\underset{\&OverBar;}{K}2}^c& \&CenterDot;\&CenterDot;& V_{\underset{\&OverBar;}{K}N}^c\end{array}\right)---\left(7\right)$

其中，N×I阶矩阵β为相关噪声场源极化参数，仅与相关噪声场源个数及观测时窗有关，β_ni(n＝1,2,...,N；i＝1,2,...,I)为其元素，表示第n个相关噪声场源第i个时窗的极化参数；K×N阶矩阵V^c为Y对应于相关噪声场源极化参数β的空间模数，仅与相关噪声场源个数及空间各测道有关，为其元素，表示第k个测道对应于第n个相关噪声场源的空间模数；ε^Y为Y中的不相关噪声；

通过对Y进行矩阵分解获得β的估计值；

第三步，利用已获取的α、β估计值以及测点信号时空数据矩阵X，求解测点对应于天然场源的空间模数U和对应于相关噪声场源的空间模数V：

测点观测数据由天然场响应和相关噪声场响应两部分叠加组成，将X写为空间矩阵与时间矩阵的乘积形式，

$X=U\mathrm{\&alpha;}+V\mathrm{\&beta;}+\mathrm{\&epsiv;}=\left[\begin{array}{cc}U& V\end{array}\right]\left[\begin{array}{c}\mathrm{\&alpha;}\\ \mathrm{\&beta;}\end{array}\right]+\mathrm{\&epsiv;}=W\mathrm{\&gamma;}+\mathrm{\&epsiv;},---\left(8\right)$

其中，W＝[UV]为总空间模数，γ＝[αβ]^*为总场源极化参数，ε为X中的不相关噪声；γ表示为：

$\mathrm{\&gamma;}={\left[\begin{array}{cc}\mathrm{\&alpha;}& \mathrm{\&beta;}\end{array}\right]}^{*}=\left(\begin{array}{cccc}{\mathrm{\&alpha;}}_{11}& {\mathrm{\&alpha;}}_{12}& \&CenterDot;\&CenterDot;& {\mathrm{\&alpha;}}_{1I}\\ {\mathrm{\&alpha;}}_{21}& {\mathrm{\&alpha;}}_{22}& \&CenterDot;\&CenterDot;& {\mathrm{\&alpha;}}_{2I}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ {\mathrm{\&alpha;}}_{M1}& {\mathrm{\&alpha;}}_{M2}& \&CenterDot;\&CenterDot;& {\mathrm{\&alpha;}}_{MI}\\ {\mathrm{\&beta;}}_{11}& {\mathrm{\&beta;}}_{12}& \&CenterDot;\&CenterDot;& {\mathrm{\&beta;}}_{1I}\\ {\mathrm{\&beta;}}_{21}& {\mathrm{\&beta;}}_{22}& \&CenterDot;\&CenterDot;& {\mathrm{\&beta;}}_{2I}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ {\mathrm{\&beta;}}_{N1}& {\mathrm{\&beta;}}_{N2}& \&CenterDot;\&CenterDot;& {\mathrm{\&beta;}}_{NI}\end{array}\right)---\left(9\right)$

由(8)式求得W以及U、V的估计值为：

[UV]＝W(11)

上角标表示共轭转置矩阵，上角标^-1表示矩阵的逆；W的展开式为：

$W=\left[\begin{array}{cc}U& V\end{array}\right]=\left(\begin{array}{cccccccc}{U}_{11}& U_{12}& \&CenterDot;\&CenterDot;& U_{1M}& V_{11}& V_{12}& \&CenterDot;\&CenterDot;& V_{1N}\\ U_{21}& U_{22}& \&CenterDot;\&CenterDot;& U_{2M}& V_{21}& V_{22}& \&CenterDot;\&CenterDot;& V_{2N}\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;& \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\&CenterDot;& \&CenterDot;\\ U_{K1}& U_{K2}& \&CenterDot;\&CenterDot;& U_{KM}& V_{K1}& V_{K2}& \&CenterDot;\&CenterDot;& V_{KN}\end{array}\right).---\left(12\right)$

$U=\left(\begin{array}{cccc}{U}_{11}& U_{12}& \&CenterDot;\&CenterDot;& U_{1M}\\ U_{21}& U_{22}& \&CenterDot;\&CenterDot;& U_{2M}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ U_{K1}& U_{K2}& \&CenterDot;\&CenterDot;& U_{KM}\end{array}\right),V=\left(\begin{array}{cccc}{V}_{11}& V_{12}& \&CenterDot;\&CenterDot;& V_{1N}\\ V_{21}& V_{22}& \&CenterDot;\&CenterDot;& V_{2N}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ V_{K1}& V_{K2}& \&CenterDot;\&CenterDot;& V_{KN}\end{array}\right)$

其中，U_km(k＝1,2,...,K；m＝1,2,...,M)为U的元素，表示第k个测道对应于第m个天然场源的空间模数；V_kn(k＝1,2,...,K；n＝1,2,...,N)为V的元素；

对第j(j＝1,2,...,J)个测点，设其各测道对应的标号分别为k+1,k+2,...,k+κ_j(k＝κ₁+...+κ_j-1)，则在空间模数矩阵U和V中，与第j(j＝1,2,...,J)个测点对应的空间模数U_(j)、V_(j)为：

$U_{\left(j\right)}=\left(\begin{array}{cccc}{U}_{(k+1)1}& U_{(k+1)2}& \&CenterDot;\&CenterDot;& U_{(k+1)M}\\ U_{(k+2)1}& U_{(k+2)2}& \&CenterDot;\&CenterDot;& U_{(k+2)M}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ U_{(k+{\mathrm{\&kappa;}}_j)1}& U_{(k+{\mathrm{\&kappa;}}_j)2}& \&CenterDot;\&CenterDot;& U_{(k+{\mathrm{\&kappa;}}_j)M}\end{array}\right),V_{\left(j\right)}=\left(\begin{array}{cccc}{V}_{(k+1)1}& V_{(k+1)2}& \&CenterDot;\&CenterDot;& V_{(k+1)N}\\ V_{(k+2)1}& V_{(k+2)2}& \&CenterDot;\&CenterDot;& V_{(k+2)N}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ V_{(k+{\mathrm{\&kappa;}}_j)1}& V_{(k+{\mathrm{\&kappa;}}_j)2}& \&CenterDot;\&CenterDot;& V_{(k+{\mathrm{\&kappa;}}_j)N}\end{array}\right),---\left(13\right)$

最后，利用求得的空间模数U和V计算各个测点的天然场阻抗张量与相关噪声场阻抗张量；对第j(j＝1,2,...,J)个测点，记天然场阻抗张量为相关噪声场阻抗张量为则

其中，分别为空间模数U_(j)、V_(j)中对应于水平磁场测道(输入端)的部分，分别为空间模数U_(j)、V_(j)中对应于水平电场测道及竖直磁场测道(输出端)的部分；具体的，对于(1)式所示的测点信号时空数据矩阵，

$U_{\left(j\right)}^H=\left(\begin{array}{cccc}{U}_{(k+1)1}& U_{(k+1)2}& \&CenterDot;\&CenterDot;& U_{(k+1)M}\\ U_{(k+2)1}& U_{(k+2)2}& \&CenterDot;\&CenterDot;& U_{(k+2)M}\end{array}\right),V_{\left(j\right)}^H=\left(\begin{array}{cccc}{V}_{(k+1)1}& V_{(k+1)2}& \&CenterDot;\&CenterDot;& V_{(k+1)N}\\ V_{(k+2)1}& V_{(k+2)2}& \&CenterDot;\&CenterDot;& V_{(k+2)N}\end{array}\right),---\left(16\right)$

$U_{\left(j\right)}^E=\left(\begin{array}{cccc}{U}_{(k+3)1}& U_{(k+3)2}& \&CenterDot;\&CenterDot;& U_{(k+3)M}\\ U_{(k+4)1}& U_{(k+4)2}& \&CenterDot;\&CenterDot;& U_{(k+4)M}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ U_{(k+{\mathrm{\&kappa;}}_j)1}& U_{(k+{\mathrm{\&kappa;}}_j)2}& \&CenterDot;\&CenterDot;& U_{(k+{\mathrm{\&kappa;}}_j)M}\end{array}\right),V_{\left(j\right)}^E=\left(\begin{array}{cccc}{V}_{(k+3)1}& U_{(k+3)2}& \&CenterDot;\&CenterDot;& V_{(k+3)N}\\ V_{(k+4)1}& U_{(k+4)2}& \&CenterDot;\&CenterDot;& V_{(k+4)N}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ V_{(k+{\mathrm{\&kappa;}}_j)1}& V_{(k+{\mathrm{\&kappa;}}_j)2}& \&CenterDot;\&CenterDot;& V_{(k+{\mathrm{\&kappa;}}_j)N}\end{array}\right).---\left(17\right)$

2.根据权利要求1所述的时空阵列差分电磁勘探方法，其特征在于，所述步骤3的第一部分中，测点信号时空数据矩阵X的构建方法为：

$X=\left(\begin{array}{cccc}{X}_{11}& X_{12}& \&CenterDot;\&CenterDot;& X_{1I}\\ X_{21}& X_{22}& \&CenterDot;\&CenterDot;& X_{2I}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ X_{K1}& X_{K2}& \&CenterDot;\&CenterDot;& X_{KI}\end{array}\right)=\left(\begin{array}{cccc}{h}_{11}& h_{12}& \&CenterDot;\&CenterDot;& h_{1I}\\ e_{11}& e_{12}& \&CenterDot;\&CenterDot;& e_{1I}\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ h_{J1}& h_{J2}& \&CenterDot;\&CenterDot;& h_{JI}\\ e_{J1}& e_{J2}& \&CenterDot;\&CenterDot;& e_{JI}\end{array}\right)---\left(18\right)$

其中，

$h_{ji}=\left(\begin{array}{c}{\left(H_{ji}\right)}_x\\ {\left(H_{ji}\right)}_y\end{array}\right),$

h_ji为第j(j＝1,2,...,J)个测点第i(i＝1,2,...,I)个时窗的水平磁场观测矢量，(H_ji)_x为其水平x方向分量，(H_ji)_y为其水平y方向分量；

当观测竖直磁场测道时，

$e_{ji}=\left(\begin{array}{c}{\left(E_{ji}\right)}_x\\ {\left(E_{ji}\right)}_y\\ {\left(H_{ji}\right)}_z\end{array}\right)$

e_ji为第j个测点第i个时窗的输出端观测矢量，(E_ji)_x为其水平x方向电场分量，(E_ji)_y为其水平y方向分量，(H_ji)_z为其竖直磁场分量；

当不观测竖直磁场测道时，

$e_{ji}=\left(\begin{array}{c}{\left(E_{ji}\right)}_x\\ {\left(E_{ji}\right)}_y\end{array}\right).$

3.根据权利要求1所述的时空阵列差分电磁勘探方法，其特征在于，所述步骤3的第二部分中，天然电磁场信号时空数据矩阵X^r的构建方法为：

当远参考点处仅布设并记录2道相互垂直的水平磁场测道时，

$X^r=\left(\begin{array}{cccc}{X}_{11}^r& X_{12}^r& \&CenterDot;\&CenterDot;& X_{1I}^r\\ X_{21}^r& X_{22}^r& \&CenterDot;\&CenterDot;& X_{2I}^r\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ X_{K^{\&prime;}1}^r& X_{K^{\&prime;}2}^r& \&CenterDot;\&CenterDot;& X_{K^{\&prime;}I}^r\end{array}\right)=\left(\begin{array}{cccc}{h}_{11}^r& h_{12}^r& \&CenterDot;\&CenterDot;& h_{1I}^r\\ h_{21}^r& h_{22}^r& \&CenterDot;\&CenterDot;& h_{2I}^r\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ h_{J^{\&prime;}1}^r& h_{J^{\&prime;}2}^r& \&CenterDot;\&CenterDot;& h_{J^{\&prime;}I}^r\end{array}\right),---\left(19\right)$

$h_{ji}^r=\left(\begin{array}{c}{\left(H_{ji}^r\right)}_x\\ {\left(H_{ji}^r\right)}_y\end{array}\right),$

为第j(j＝1,2,...,J')个远参考点第i(i＝1,2,...,I)个时窗的水平磁场观测矢量，为其水平x方向分量，为其水平y方向分量；

当远参考点处增加布设并记录2道相互垂直的水平电场测道时，

$X^r=\left(\begin{array}{cccc}{X}_{11}^r& X_{12}^r& \&CenterDot;\&CenterDot;& X_{1I}^r\\ X_{21}^r& X_{22}^r& \&CenterDot;\&CenterDot;& X_{2I}^r\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ X_{K^{\&prime;}1}^r& X_{K^{\&prime;}2}^r& \&CenterDot;\&CenterDot;& X_{K^{\&prime;}I}^r\end{array}\right)=\left(\begin{array}{cccc}{h}_{11}^r& h_{12}^r& \&CenterDot;\&CenterDot;& h_{1I}^r\\ e_{11}^r& e_{12}^r& \&CenterDot;\&CenterDot;& e_{1I}^r\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ h_{J^{\&prime;}1}^r& h_{J^{\&prime;}2}^r& \&CenterDot;\&CenterDot;& h_{J^{\&prime;}I}^r\\ e_{J^{\&prime;}1}^r& e_{J^{\&prime;}2}^r& \&CenterDot;\&CenterDot;& e_{J^{\&prime;}I}^r\end{array}\right)---\left(20\right)$

其中

$e_{ji}^r=\left(\begin{array}{c}{\left(E_{ji}^r\right)}_x\\ {\left(E_{ji}^r\right)}_y\end{array}\right)$

为第j(j＝1,2,...,J')个远参考点第i(i＝1,2,...,I)个时窗的水平电场观测矢量，为其水平x方向分量，为其水平y方向分量。

4.根据权利要求1所述的时空阵列差分电磁勘探方法，其特征在于，所述步骤3的第三部分中，设第k(k＝1,2,...,K)个差分磁场测道由第j和j'(j,j'＝1,2,...,J且j≠j')两个测点的d方向(d方向为x、y两个方向之一)水平磁场进行差分所得，则第k个差分磁场测道第i个时窗的水平磁场差分数据Y_ki为：

Y_ki＝H_j,d,i-H_j',d,i,(21)

其中，H_j,d,i、H_j',d,i分别为第j和j'个测点的d方向第i个时窗的水平磁场观测数据。

5.根据权利要求1所述的时空阵列差分电磁勘探方法，其特征在于，所述步骤4的第一步和第二步中，矩阵分解方法包括奇异值分解、主成分分析法和稳健性主成分分析法。

6.根据权利要求1所述的时空阵列差分电磁勘探方法，其特征在于，所述步骤4的第一步中，采用的矩阵分解方法为奇异值分解；具体步骤为：

[U,S,V]＝svd(X^r),(22)

U^r＝(U)^K'×M,α＝(SV ^*)^M×I,(23)

其中，svd表示矩阵的奇异值分解，U,S,V分别为奇异值分解所得的参数矩阵，上角标^*表示矩阵转置；(U)^K'×M表示U的1～M列组成的K'×M阶矩阵，(SV ^*)^M×I表示SV ^*的1～M行组成的M×I阶矩阵。

7.根据权利要求1所述的时空阵列差分电磁勘探方法，其特征在于，所述步骤4的第二步中，采用的矩阵分解方法为奇异值分解；具体步骤为：

[U,S,V]＝svd(Y),(24)

V^c＝(U) ^K×N,β＝(SV ^*)^N×I,(25)

其中，svd表示矩阵的奇异值分解，U,S,V分别为奇异值分解所得的参数矩阵，上角标^*表示矩阵转置；(U) ^K×N表示U的1～N列组成的K×N阶矩阵，(SV ^*)^N×I表示SV ^*的1～N行组成的N×I阶矩阵。