Home Journals Progress in Geophysics
Progress in Geophysics

Abbreviation (ISO4): Prog Geophy      Editor in chief:

About  /  Aim & scope  /  Editorial board  /  Indexed  /  Contact  / 

Anti-interference processing for CSAMT based on deep learning and joint de-noising

  • WeiQiang LIU , 1 ,
  • PinRong LIN , 2, * ,
  • RuJun CHEN 3 ,
  • Kun ZHANG 1 ,
  • ChangXin CHEN 1 ,
  • Xu LIU 4
Expand
  • 1 SinoProbe Laboratory, Chinese Academy of Geological Sciences, Beijing 100037, China
  • 2 Key Laboratory of Geophysical Electromagnetic Probing Technologies of Ministry of Natural Resources, Institute of Geophysical and Geochemical Exploration of Chinese Academy of Geological Science, Langfang 065000, China
  • 3 Institute of Geosciences and Info-Physics, Central South University, Changsha 410083, China
  • 4 School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China

Received date: 2023-09-03

  Online published: 2024-09-29

Copyright

Copyright ©2024 Progress in Geophysics. All rights reserved.

Abstract

Controlled-Source Audio Magnetotelluric (CSAMT) is a near-surface geophysical method that developed on the basis of Magnetotelluric method (MT). With the development of social economy,the data quality of CSAMT has also been seriously disturbed by noise interference. In practical exploration,the time series of electromagnetic field is usually superimposed with large-scale trend drift,short-term sudden strong interference and peak impulsive outliers,resulting in the distortion of the calculated resistivity spectrum. In this paper,an anti-interference processing method based on deep learning and joint de-noising is proposed to preprocess CSAMT time series. Firstly,a forward algorithm of electromagnetic time series of layered earth controllable source is proposed,which is used to generate standard electromagnetic signals without noise interference. Then,a Long and Short Term Memory Neural Network (LSTM) classifier is trained to recognize the noise. Finally,the improved Empirical Mode Decomposition (EMD) algorithm,correlation based data selection algorithm and robust statistical algorithm are jointly used to de-noise the CSAMT time series. The test results by simulated data show that the recognition accuracy of LSTM for noise interference can reach more than 95%,and the three noise reduction algorithms can reduce the data error from about 20% to less than 3%. Finally,the proposed method is applied to the actual data set of a metal mining area in Inner Mongolia. the accuracy of low-frequency resistivity and phase was effectively improved.

Cite this article

WeiQiang LIU , PinRong LIN , RuJun CHEN , Kun ZHANG , ChangXin CHEN , Xu LIU . Anti-interference processing for CSAMT based on deep learning and joint de-noising[J]. Progress in Geophysics, 2024 , 39(4) : 1457 -1473 . DOI: 10.6038/pg2024HH0341

0 引言

可控源音频大地电磁法(CSAMT)是在天然源大地电磁法(MT)的基础上,为了进一步增强信噪比而提出的一种以人工发射大功率电磁波来代替天然电磁场源的频率域电磁勘探方法.早期主要是模仿MT的工作方式,在电磁波传播的远区(平面波区)进行数据采集,采用与MT类似的手段进行数据处理和反演解释.近年来,随着仪器装备和计算方法的进步,可控源音频大地电磁法已经可以在全区进行数据采集、处理和解释(汤井田和何继善,2005底青云和王若,2008何继善,2010).虽然可控源音频大地电磁法已经在国内外都得到了普及,数据采集和反演解释技术得到了长足发展,但随着社会经济的迅速发展,各类电磁干扰日益严重,影响了其数据质量.噪声干扰的压制需要从可控源时间序列入手,分析噪声干扰特点,根据信号特征,研发抗干扰数据处理技术.
目前在电磁勘探中,天然源大地电磁法的信号处理与抗干扰技术最为成熟,除了常规的频谱分析、阻抗估计、数据筛选等方法外,国内外学者先后将数学形态滤波、稀疏分解、小波分析、子空间增强、压缩重构、递归分析等新算法引入大地电磁处理中,提高了数据质量.相比于大地电磁法,很长一段时间内,可控源音频大地电磁法主要通过增大发射电流功率和改进测量方式来提高信号质量, 目前国内自主研制的可控源发射机功率已达到200 kW以上(何展翔等,2020Zhou et al., 2018王猛等,2019).随着硬件本身的限制,进一步增大发射功率的难度越来越大.而且随着发射机功率的增大,能源消耗会进一步增大、施工安全性进一步降低.所以开展抗干扰数据方法研究依然十分有必要.
近年来,国内外学者关于可控源音频大地电磁法抗干扰处理的研究主要有以下几类:
首先是频率域的降噪算法.在天然源大地电磁中应用较广的频谱挑选、稳健叠加、远参考等算法也可应用于可控源音频大地电磁数据处理.李予国和段双敏(2014)将傅里叶变换、预白处理、漂移矫正、多次叠加等技术用于海洋可控源数据的预处理.刘宁(2015)提出通过频率域振幅叠加、时变双边滤波等方法压制低频的天然电磁场和海底洋流噪声.杨洋等(2018)在频域对主频及各谐波成分进行Hilbert分解,根据包络线来估计信噪比来筛选数据.Li等(2021a)进一步提出将快速傅里叶变换、互补集合模态分解和移不变稀疏编码相结合用于压制工频干扰、基线漂移、环境噪声等人文干扰.Jiang等(2022)将引导图像滤波用于可控源音频大地电磁视参数的静态效应校正和噪声干扰压制.频率域处理方法主要是基于傅里叶变换进行频谱挑选,当噪声干扰与电磁信号的频谱混叠严重时,应对能力有限.
其次是时间域的数据处理算法.学者们先后提出了自适应双向均方差阈值法、灰色系统识别方法、稳健中值滤波、有理函数滤波等改进算法对电磁勘探时间序列中的粗大误差进行判别和剔除(雷达等,2010Myer et al., 2011张必明等,2015Ge et al., 2016Mo et al., 2017陈超健等,2019崔扬等,2023胡艳芳等,2024).此外,Jiang等(2018)采用柯西振荡蛙跳算法优化的小波包对视电阻率数据进行降噪以提高反演效果.Li等(2021b)提出以归一化互相关系数、振幅绝对最大值和趋势波动分析为指标,通过模糊C均值聚类对时间序列中的强干扰数据段进行识别和剔除,提高了数据处理的自动化水平.Zhang等(2022)进一步提出通过灰狼优化算法(IGWO)和支持向量机(SVM)对强干扰数据段进行识别,并利用数字相干技术提取有效频点的频谱幅度,从而重构原始信号.目前已有的时域降噪算法对于粗大误差和间断性强干扰效果较好,对全时段持续出现的多源混合噪声降噪能力依然有限.
还有一种降噪思路是采用等效源方法(MacLennan and Li, 2013; Tang et al., 2018),把噪声干扰、静态效应、电流畸变等干扰当作信号源的一种进行处理和分析.周聪等(2019)进一步提出了时空阵列混场源电磁法理论,通过求解系统响应将不同场源的信号分离以获取高质量的大地电磁和可控源音频大地电磁信号.上述思路具有较好的降噪前景但目前尚不成熟,需要进一步开展实测数据测试.
近年来,以人工神经网络为代表的机器学习、深度学习等智能算法得到迅速发展,在大地电磁信噪识别、环境背景噪声分类、瞬变电磁数据降噪、天然地震信号识别等领域得到初步应用,前人的研究表明,深度学习方法对在信噪识别和分离方面具有更高的精度和效率(Meier et al., 2019; Cao et al., 2019; Wu et al., 2022; Li et al., 2019, 2023).
前人的研究有效提高了可控源电磁勘探数据质量,但已有的信号处理手段主要还是基于电磁场频谱幅值进行数据挑选或多次叠加,而且多是针对噪声干扰统计特征明确的数据.在实际勘探中,各类电磁干扰种类多样、变化复杂且持续时间较长,靠单一方法很难挑选出高品质的数据或频谱.为进一步提高可控源音频大地电磁数据质量,本文着重从以下两点开始研究,一是深入分析可控源信号及各类电磁干扰的时间序列,时频域相结合,通过深度学习识别噪声类型;二是对多种有针对性的信号处理手段进行优化改进构建联合降噪方法,在识别各种噪声干扰后,自动选择相应的信号处理技术进行压制,提高降噪效率和精度.

1 基于深度神经网络的CSAMT信噪识别

1.1 电磁信号与噪声干扰样本生成

获得不含噪声干扰的参考信号是开展信噪识别和分离的前提,由于获取大量高质量电磁信号的成本很高,本文提出了一种一维全波形正演算法来生成CSAMT电磁信号作为参考样本.人工源电磁勘探中常用的发射电流波形包括矩形波、双极性波、双频波、2n序列伪随机波、m序列扩频波等,都属于周期性信号(Ziolkowski et al., 2011; Zhao et al., 2013; Song et al., 2018; Hu et al., 2020).根据傅里叶变换理论,周期信号可展开成正余弦级数求和的形式.通过假设地下介质为均匀半空间或层状介质,并将发射电流分解为不同频率的正余弦波级数之和,正演出单频电磁响应后再组合叠加可得到与实测信号类似的时间序列波形.
周期性电流信号的傅里叶级数分解公式为(Zahn and Roskies, 1972):
其中各参数的表达式如下:
在确定发射电流和采集布置后,通过数值滤波或积分变换(Key, 2012)可以模拟层状大地在不同频率下的表面电磁场响应,公式简化如下:
其中,F为感应电场或磁场,h为层厚度,ρ为每层的电阻率,ω为发射电流的角频率.计算出单频的电磁响应之后,通过傅里叶级数叠加将所有的单频响应相加,便可以获得时域中的电磁时间信号:
原则上需要无穷多组频率响应的叠加才能得到真实的时间序列.但由于发射电流的谐波能量随着频率的增加而急剧衰减,因此公式(1)和公式(4)中的无穷组分量可以通过数千甚至数百个分量的叠加来近似.通过使用并行计算,可以快速上万组与测量数据相似的标准电磁场序列.本文采用1~10000组傅里叶级数来合成上述波形,并计算了合成波形与标准波形的均方相对误差,当傅里叶级数超过5000项时,方波和双极性波拟合误差已低于0.1%, 2n序列伪随机波和m序列扩频波拟合误差已低于1%.图 1展示了不同类型电流波形和频谱及拟合误差随傅里叶级数数量的变化曲线.
图1 人工源电磁勘探常用激励电流波形和频谱及拟合误差随傅里叶级数数量的变化曲线

从上至下分别为方波、双极性波、扩频波、伪随机波等.

Fig 1 Excitation current waveform and its frequency spectrum commonly used in artificial source electromagnetic exploration, the changing of waveform reconstruction error with the increase of the number of Fourier series is also shown

From top to bottom are square wave, bipolar wave, spread spectrum wave, pseudo-random wave, etc.

本文主要研究陆地CSAMT勘探中常见的三种噪声干扰,包括低频趋势项漂移、短时突发强干扰和尖峰脉冲干扰.低频趋势项漂移由大地自然电流、工业杂散电流、电极极化电位和仪器零点漂移等引起,对许多电磁勘探方法的低频数据都有很大的影响.短时突发强干扰由突发性的工业设施和人类活动引起,其幅值远大于正常信号,且变化剧烈.尖峰脉冲干扰由天电干扰、大功率机械开关等引起的,是时域内的振荡峰值.在实际勘探中,上述三种干扰可能单独出现,也可能耦合出现,可能只出现在某个时间段中,也有可能叠加在整个时间序列上.电场和磁场都会受到上述噪声干扰的影响,导致计算的视电阻率和视相位发生畸变.只有预先抑制时间序列中的噪声,才能得到准确的阻抗估计结果.本研究通过数学模拟生成三种干扰的训练样本.通过正弦函数、多次函数、线性函数和指数函数的叠加模拟低频趋势干扰.通过间断出现的强高斯信号模拟短时突发强干扰,该噪声的整体振幅大于正常信号,但仅在某些时间段连续出现.尖峰脉冲干扰由振幅远大于正常信号且在整个周期内零散出现的离群值模拟.然后,将上述干扰随机组合并添加到纯电磁信号中,三类干扰随机组合共有八种类型,对应不同的降噪处理方法.

1.2 长短时记忆神经网络分类器

本文通过长短时记忆神经网络(LSTM)对电磁时间序列进行分类,进而识别噪声干扰类型.长短时记忆神经网络(LSTM)由常规的误差后向传播神经网络(BPNN)和循环神经网络(RNN)发展而来,是时间序列处理中常用的一种深度学习方法(许滔滔等,2020).长短时记忆神经网络(LSTM)利用不同长度的记忆模块代替普通的隐含层节点,在传递不同时间步信息的同时,可有效避免常规神经网络在输入样本时间步过多时出现的梯度爆炸、不易收敛等问题.从而利用较小的网络规模高效地学习到长序列数据信息.
与常规神经网络类似,LSTM同样包括输入层、隐含层、输出层.LSTM神经网络输入样本是时间序列不同时间步的统计特征.输入层对原始输入进行预处理,然后通过归一化将其转换为与隐藏层相匹配的标准化输入数据.隐含层由许多LSTM网络节点组成,用于构造从输入到输出的非线性映射.输出层包括全连接层、SOFTMAX层和分类层,它们共同将隐藏层的输出转换为时间序列的类别标签.上述三个网络层都由一系列要优化的权重系数组成:
本文使用仿真生成的含噪声干扰的电磁时间序列作为输入,噪声干扰类型标签为输出,对上述神经网络系统进行训练,以获得最佳权重系数.在训练中,预测输出可以简化表示为:
其中F是一个显式非线性变换函数,也称为激活函数.通过最小化预测的输出YP和实际的输出YT之间的差异φ,可以获得最佳权重函数:
在迭代过程中,Adam优化器(Khobragade et al., 2020)用于优化每个迭代过程中的参数,直到目标函数趋于收敛.
与传统神经网络所不同的是,LSTM在每个隐含层的节点上设置了两个传输参数.这样既可以考虑不同时间步长的特点,又可以避免数据积累引起的梯度爆炸.本文设计的网络系统共包括两个LSTM隐含层,每层50个网络节点.图 2展示了LSTM每个网络节点的结构.该节点包括两个传输单元和.当前时间步的传输单元ht, Ct受当前输入和前一时间步的传输单元ht-1, Ct-1的影响.计算公式如下:
图2 一个LSTM节点的结构,当前时间步的传输单元ht, Ct受当前输入xt和前一时间步长传输单元ht-1, Ct-1的综合影响

(a)原始电场和拟合电场;(b)原始磁场和拟合磁场;(c)含噪电场和拟合电场;(d)含噪磁场和拟合磁场.

Fig 2 The structure of an LSTM node, the transmission units ht and Ct of the current time step is affected by the current input xt and the transmission units ht and Ct of the previous time step

其中,ft是一个遗忘门,用于确定前一个单元状态Ct-1保留到当前时间步的比例;it是一个输入门,用于确定新增量ΔCt保留到当前时间的比例;ot是一个输出门,用于确定单元状态Ct作为当前输出值的比例.所有参数ot,ΔCtitft均根据ht-1xt进行计算.公式如下:
其中,sigmoid()和tanh()为激活函数,用于将输入变量进行非线性变换.在生成时间步的传输单元状态和增量时,采用tanh()函数,其输出在-1~1之间,符合大多数场景下传输单元特征分布以0为中心的特点,且tanh()函数在输入为0附近时具有更大的梯度,可使模型快速收敛;生成输入门、输出门和遗忘门时,采用sigmoid()函数,因为其输出在0~1之同,符合门控的物理定义,且当输入较大或较小时,其输出会非常接近1或0,从而保证该门开或关.在上述公式中,Wi, Wf, WC, 和Wo是需要估计的节点权重.
本文以电磁时间序列全时段的时频域统计特征为输入,是否含有上述三种噪声干扰的不同组合为噪声类型标签,训练神经网络分类器.当训练样本的整体拟合误差低于误差容限时,训练结束.最后将训练得到的网络系统用于实际信号噪声识别.为增强LSTM网络的稳定性,采取了以下改进措施:
(1) 提取时间序列的统计特征作为输入样本.在分类问题中,如果直接以完整时间序列为输入样本,则所需的神经元节点数及样本数将非常巨大.本研究先对含噪声干扰的时间序列进行分段,并统计得到不同时间段内表征信号类型的统计特征.时域主要是通过统计分析方法,对电磁时间序列在不同时间步的均值、标准差、相关度、离群值比例、曲线光滑度等统计特征进行汇总.频域主要是通过时频分析算法进行统计,统计原始数据在不同时间段的基频振幅、基频相位、主频能量比例等特征.对于不含
噪声干扰的纯电磁信号与含有不同类型噪声干扰的混合信号,其统计特征在时域或频域必然存在不同点.最后以上述时频域统计特征为判别指标,通过长短时记忆神经网络(LSTM),对电磁信号和噪声干扰进行分类识别.
(2) 通过交叉验证、随机扰动、并行计算等手段增强网络的稳定性.深度神经网络由于网络层数多、结构复杂,容易出现模型稳定性差、易陷入过拟合等问题.本文通过施加随机扰动、交叉验证、并行计算等方式提高神经网络训练效果.即同时对所有原始输入样本加入随机扰动,然后训练,以提高网络稳定性.同时在模型训练过程中就将输入的训练样本随机分为训练、测试和验证三部分,用于训练、验证和测试的样本比例分别为80%、10%和10%.在每次训练的同时就进行测试,既要求模型的数据拟合差最小,又要求模型能充分适应新数据.由于交叉验证是在训练过程中随机选取数据,所以相同数据多次训练所得的模型略有不同.本文采用并行技术同时训练多个模型,从中选择数据测试误差最小的神经网络作为最终结果来处理新数据.

2 基于联合降噪算法的CSAMT信噪分离

由于可控源音频大地电磁时间序列中含有多种类型的噪声干扰,单一方法难以有效压制混合噪声,本文采用联合降噪算法对CSAMT时间序列进行处理,包括改进经验模态分解算法、相关分析数据挑选算法和稳健统计叠加算法(如图 3所示).上述三种方法是在传统方法的基础上,结合可控源电磁信号特点进行优化改进得到.首先,传统的经验模态分解算法可用于周期性和非周期性信号的去趋势项,但计算过程需要两层循环,计算效率较慢,考虑到可控源电磁信号为周期性信号,本文对传统经验模态分解算法进行改进,保证去趋势项精度的同时,将两层循环变为一层循环,降低了计算成本.同时,传统的相关分析数据挑选算法只用于评价和剔除不同信号段数据质量,未对强干扰数据段做进一步处理,本文将可控源电磁时间序列正演和相关分析相结合,用正演仿真的信号去拟合和代替干扰数据段,从而实现了电磁信号的恢复,减少了数据损失.最后,考虑到常规均值叠加方法仅能压制高斯随机噪声,无法应对尖峰脉冲干扰,本文采用稳健统计叠加方法对多周期可控源电磁数据进行处理,并优选了权重函数,可同时压制高斯随机噪声和尖峰脉冲干扰.
图3 CSAMT电磁场时间序列联合降噪方法流程

Fig 3 Joint de-noising process of CSAMT electromagnetic field time series

2.1 改进经验模态分解算法

为分离可控源音频大地电磁时间序列中的低频趋势项漂移,提出了一种改进的经验模态分解方法.对于周期为T的信号,一个周期内可能存在多个局部固有极大值和极小值.一个周期内不同位置的局部极值点可能不相等,但不同时期同一位置的极值点应该相等.趋势项漂移的存在将改变周期信号之间的等效性.本文以周期信号之间的周期等效性为准则,对信号进行分解,确定趋势项,具体步骤如下:
(1) 首先通过对多周期可控源音频大地电磁时间序列进行叠加,以获得一个周期的相对高质量的数据.将叠加后的数据进一步分成n段,分别计算每段的最大值和最小值.这样就得到了每一个周期内部n组局部极大值和极小值.
(2) 然后依次连接不同周期在同一位置的极值点,通过插值得到n条完整的上下包络线.由于局部极值通常不是位于时间序列的末端,而是位于时间序列的内部,因此在时间序列两端存在少量包络线无法覆盖的数据.上述问题称为端点效应.对于具有大量周期的高频段电磁数据,包络线外的数据可以直接删除,不用于后续处理.对于周期数较少的低频段数据,本文通过二次函数近似拟合前两个周期和最后两个周期的趋势项,从而对包络线无法覆盖的数据进行处理.
(3) 最后对所有的包络进行叠加,得到趋势项漂移的初步估计.然后从原始数据中减去估计的趋势项,得到新数据,并对新数据重复上述操作(1)和(2).原则上当原始时间序列包含的趋势项被完全分离后,不同周期同一位置的采样点应大致相等,即时间序列的周期性得到恢复.因此,本文计算所有周期同一位置数据点的标准偏差之和来确定收敛阈值,公式如下:
其中,N是周期数,M是一个周期的数据点数.两次迭代的标准差变化率定义如下:
然后,当ε < 1/100时,迭代停止.此时,残留的趋势项漂移近似为一条均值为0的直线,电磁信号的周期性得到恢复.
此外,为了确定步骤(1)中的参数n,本文采用模拟数据对其进行了测试.将30个周期的矩形波信号(每个周期300个数据点)作为不含趋势项干扰的标准信号,将不同频率的正弦波组合作为趋势漂移干扰,将二者叠加得到含趋势项干扰的仿真信号用于后续测试.同时,在时间序列中加入高斯随机噪声和尖峰脉冲噪声,以测试去趋势项算法的稳定性.将参数n依次设置为1~300,计算拟合的趋势.同时计算拟合的趋势项与实际趋势项的相对误差以及时间成本.仿真结果如图 4所示,随着参数n的增加,时间成本呈线性增加.同时,拟合误差随着参数n的增加波动性减小.这是因为当n较大时,可以获得更多的包络线,并且可以通过多组包络线的叠加来抑制其他噪声干扰的影响.然而,当n大于数据点数的一半时,计算误差明显反弹,最终误差没有减小.因此,将n值设置为少于数据点数量的一半是合理的.图 5显示了此次仿真的矩形波信号、正弦趋势漂移信号、包含趋势项的信号,以及参数n为50时的拟合趋势.真实趋势漂移和拟合趋势漂移几乎完全一致,表明去趋势项算法效果较好.
图4 在改进的经验模态分解算法中,拟合误差和时间成本随参数n的变化

Fig 4 Variation of fitting error and time cost with parameter n in the modified empirical mode decomposition algorithm

图5 用改进的经验模态分解算法从扰动数据中分离趋势漂移的示意图

Fig 5 Schematic diagram of separating trend drift from disturbance data with modified empirical mode decomposition algorithm

2.2 相关分析数据挑选算法

在可控源音频大地电磁探测中,每个测点的数据都经历了很长的采集时间.在某些时间段,可能会出现突发强噪声干扰,淹没电磁信号.本文进一步使用相关分析和波形拟合来处理这种干扰.
(1) 在实际勘探中,虽然采集的电磁场波形是可变的,但它们都与发射电流相关.根据周期数将实际采集的电磁场数据分为N段,每段有m个数据点.然后,分别计算每个数据段电磁场数据Uk与同步发射电流Ik之间的皮尔逊相关系数c(k)(Benesty et al., 2008):
当电磁场受到突发强噪声干扰时,与发射电流的相关性会急剧下降.因此,相关系数被用作评价数据质量的指标.
(2) 根据N个数据段的相关系数分布,将阈值设置为0.5和0.25,以选择高质量数据.当相关系数大于0.5时,数据段通常是高质量的,可以直接用于后续叠加.当相关系数小于0.25时,说明数据段受噪声干扰而严重失真,应予以消除.当相关系数在0.25和0.5之间时,使用20层地球模型计算的拟合波形来替代该数据段.
(3) 通过优化方法计算出与实测数据最接近的拟合波形.目标函数如下:
式中,To是实际观测到的时间序列,Tf是用公式(1)~(4)正演模拟得到的时间序列.h1, h2, …,hn-1ρ1, ρ2, …,ρn是分层模型各层的电阻率和厚度参数.通过迭代,可以得到最佳的模型参数和拟合波形.本文使用的优化算法是内点法(Nocedal et al., 2014).
本文进一步用一组实测的电场和磁场数据来说明波形拟合的效果.图 6ab显示了原始电磁场和拟合结果.拟合波形近似恢复实测数据段.图 6cd显示了添加强随机干扰后的信号和拟合结果.电场和磁场的信噪比(SNR)分别为-1.2和-1.7.虽然噪声干扰很强,但拟合信号仍具有较高的精度.图 7ac为原始电磁场和拟合电磁场的低频频谱图,原始频谱和拟合频谱近似相等.图 7bd分别为受干扰电磁场的低频频谱图,噪声对电场和磁场产生频谱畸变,相对误差的均方根分别为25.17%和44.96%. 然后对拟合波形的频谱进行计算,误差分别减小到3.44%和7.78%.上述仿真结果表明,通过波形拟合可以在一定程度上恢复强噪声干扰下的电磁场信号.
图6 原始电磁场和波形拟合电磁场对比示意图

Fig 6 Comparison diagram of original electromagnetic field and fitting electromagnetic field

(a) Original electric field and fitting electric field; (b) Original magnetic field and fitting magnetic field; (c) Noisy electric field and fitting electric field; (d) Noisy magnetic field and fitting magnetic field.

图7 原始电磁场、受污染电磁场和拟合电磁场的低频频谱(只有前19个奇谐波)

(a)原磁场和拟合磁场的频谱;(b)原磁场和拟合磁场的频谱;(c)原磁场、污染磁场和拟合磁场的频谱;(d)原磁场、污染磁场和拟合磁场的频谱.

Fig 7 The low-frequency spectrum (only the first 19 odd harmonics) of the original, contaminated, and fitting electromagnetic fields

(a) The spectrums of original and fitting electric fields; (b) The spectrums of original and fitting magnetic fields; (c) The spectrums of original, contaminated and fitting electric fields; (d) The spectrums of original, contaminated and fitting magnetic fields.

2.3 稳健统计叠加算法

通过改进经验模态分解消除大尺度漂移,采用相关分析数据挑选去除强干扰数据后,需要对剩余的多周期数据进行叠加,得到单周期高质量信号,进而计算视电阻率谱.本文在后续信号处理中采用稳健统计叠加方法(Huber,1996),在压制高斯随机噪声的同时去除由尖峰脉冲干扰引起的离群值.稳健统计原理如下:对于重复观测的数据点{yi}:
其中,ξ1, ξ2, ξ3, …, ξN为误差序列,N为总观测次数.稳健统计方法是根据多次观测值来估计真实值,根据最大似然估计准则,稳健均值可以通过极小化以下目标损失函数来实现:
式中,θ为“位置参数”,代表观测数据的真实值.σ是“尺度参数”,代表观测值和实际值之间的残差分布.该函数R()为似然函数.上述最优化问题可以通过求解影响函数方程来计算:
迭代算法可用于求解公式(16)中的位置参数.首先,将原始观测序列的中值和中值绝对偏差作为初始估计值:
然后进行迭代,迭代公式定义为:
其中ψ′是ψ的导数.ψ是影响函数,有多种类型,包括:“Andrews”, “Beaton”, “Talwar”, “Cauchy”, “Welsch”, “Huber”, “logistic”, “fair”, “Hampel”等(Holland and Welsch, 1977).由于可控源音频大地电磁场是多周期信号,周期数通常为几十到几百个,去除少量离群值不会影响计算结果.因此,本文选择“Talwar”权重函数.公式如下:
“Talwar”是一个硬阈值损失函数,异常值的权重为0.对于周期数较少的信号,例如低频激电信号,可以使用Huber和其他权重函数为异常值赋予较低的权重,而不是直接去除.

3 仿真数据降噪测试与误差分析

本文随机生成了1250组五层介质模型.各层厚度在50~500 m之间,各层电阻率在10-1~105 Ω ·m之间.发射电流为15组不同基频的矩形波,其基频依次为0.5 Hz、1 Hz、2 Hz、4 Hz、8 Hz、16 Hz、32 Hz、64 Hz、128 Hz、256 Hz、512 Hz、1024 Hz、2048 Hz、4096 Hz、8192 Hz.每个层状模型先正演生成电磁场频率响应,然后通过前述的傅里叶级数合成方法生成相应的电磁场时间序列.然后,在纯时间序列中加入仿真生成的八种组合噪声干扰.每个频率总共生成了20000组含噪时间序列(电场和磁场)作为训练样本.然后将上述样本按80%, 10%, 10%的比例随机分为训练集、交叉验证集、测试集,即样本数分别为16000,2000,2000.输入数据为各时间序列在不同时间步的时频域统计特征,输出数据为混合噪声干扰类型标签.训练过程中交叉验证精度为94.60%.训练完成之后进行测试,噪声类型识别精度为95.35%, 达到预期要求.
为了进一步检验噪声识别的准确性,使用一组由两层大地模型正演得到与训练样本不同的电磁时间序列作为标准信号,然后添加仿真噪声,对训练后的LSTM分类器和三种降噪方法进行测试.在该模型中,第一层的电阻率为1000 Ω ·m,层厚为200 m,第二层为均匀半空间,电阻率为100 Ω ·m.发射电流源位于坐标原点,观测点位于x=5000 m和y=5000 m处.发射电流为一组矩形波,基频为0.5 Hz、1 Hz、2 Hz、4 Hz、8 Hz、16 Hz、32 Hz、64 Hz、128 Hz、256 Hz、512 Hz、1024 Hz,2048 Hz、4096 Hz、8192 Hz等.图 8显示了三组发射电流、感应电场和感应磁场数据(基频分别为0.5 Hz、32 Hz和256 Hz).在时间序列中随机添组合的噪声干扰,共生成1000组含噪声干扰的电磁场数据.经测试,LSTM对电场和磁场的识别准确率分别为95.30%和96.23%.
图8 三个高质量电磁场的时间序列

(a)、(b)和(c)分别是基频0.5 Hz、32 Hz和256 Hz的发射电流;(d)、(e)和(f)分别是基频0.5 Hz、32 Hz和256 Hz的电场;(g)、(h)和(i)分别是基频0.5 Hz、32 Hz和256 Hz的磁场.

Fig 8 Time series of three high-quality electromagnetic fields

(a), (b) and (c) are emission currents with fundamental frequencies of 0.5 Hz, 32 Hz and 256 Hz, respectively; (d), (e) and (f) are electric fields with fundamental frequencies of 0.5 Hz, 32 Hz and 256 Hz, respectively; (g), (h) and (I) are magnetic fields with fundamental frequencies of 0.5 Hz, 32 Hz and 256 Hz, respectively.

本文进一步展示了三种降噪算法的处理流程.计算了基频和三次谐波频率下的视电阻率和视相位.图 9显示了添加模拟噪声干扰(基频32 Hz)后电磁场的波形.由于空间限制,未显示其他14个频率点的波形.
图9 对两层大地模型正演的电磁场添加噪声干扰后的波形(基频为32 Hz)

(a)发射电流;(b)感应电场;(c)感应磁场.

Fig 9 Waveform after adding noise interference to the electromagnetic field forward simulated by the two-layer geodetic model (the fundamental frequency is 32 Hz)

(a) Emission current; (b) Induced electric field; (c) Induced magnetic field.

然后,采用三种方法对15组时间序列进行处理,得到了30个频点的视电阻率和视相位曲线.第一种方法是通过均值叠加计算多周期电场和磁场的平均值,然后计算视参数.第二种方法是仅用稳健统计方法对电场和磁场的多周期时间序列进行叠加,然后计算视参数.第三种方法是先用改进经验模态分解和相关分析方法分离大尺度趋势漂移和短时强干扰,然后用稳健统计方法叠加剩余数据,计算视参数.将三种方法的结果与无噪声干扰的数据进行了比较,如图 10所示.从图 10ab可以看出,噪声的存在使均值叠加的计算结果失真,曲线上有许多跳点.视电阻率和视相位与实际值的均方相对误差分别为26.58%和15.02%.从图 10cd可以看出,仅通过使用稳健统计方法进行处理的数据失真没有减少,而是进一步增加.视电阻率和视相位的相对均方误差分别为26.11%和16.37%.这主要是因为原始时间序列不仅包含高斯随机噪声和尖峰脉冲噪声,还包含大尺度趋势漂移和短时突发强干扰.从图 10ef可以看出,噪声干扰得到了有效抑制,视电阻率和视相位曲线得到了平滑.计算的均方相对误差分别为1.29%和0.63%, 数据质量显著提高.
图10 通过三种方法得到的视电阻率和视相位,并与真实值作对比

(a)和(b)为均值叠加法获得的视电阻率和相位;(c)和(d)为稳健叠加法获得的视电阻率和相位;(e)和(f)为联合降噪法获得的视电阻率和相位.

Fig 10 Apparent resistivity and apparent phase obtained by three methods and compared with the real values

(a) and (b) are the apparent resistivity and phase obtained by the mean superposition method; (c) and (d) are the apparent resistivity and phase obtained by the robust superposition method; (e) and (f) are the apparent resistivity and phase obtained by the joint noise reduction method.

4 某矿区实测电磁数据测试与分析

为了进一步对联合降噪方法进行验证,在某测区实测数据中挑选单周期高质量的电场、磁场及发射电流数据作为不含噪声干扰的标准信号进行测试.首先将原数据扩展到50~200个周期,然后依次仿真生成低频趋势项干扰、短时突发强干扰、尖峰脉冲干扰、高斯随机噪声等,再通过随机组合叠加产生多源混合噪声.然后将仿真得到的噪声干扰添加到高品质电场和磁场信号中,进行信噪识别和分离.
图 11展示了基频32 Hz电场和磁场时间序列添加仿真噪声干扰后的波形,受篇幅限制,其余14个频点的波形没有展示.然后仍然采用与前述仿真数据处理类似的三种方法对15组时间序列进行处理,分别是均值叠加方法、仅采用稳健统计叠加方法和采用联合降噪方法.三种方法所得结果分别与不含噪声干扰数据的计算结果进行了对比.由于实测数据的最高采样率为24000 Hz,低于最后两个方波的三次谐波频率的两倍,不满足采样定律(Por et al., 2019),因此本文仅计算了前28个频率点的视电阻率和视相位曲线.最终计算结果如图 12所示.由图 12ab可见,噪声干扰的存在使得均值叠加计算结果产生了畸变,曲线上产生很多跳点,计算得到视电阻率和视相位与真实值的均方相对误差分别为17.06%与8.75%.由图 12cd可见,仅采用稳健统计方法进行处理,数据畸变不仅没有减小,反而进一步增大,计算视电阻率和视相位均方相对误差分别为19.26%和19.39%, 这主要是因为原始时间序列中不仅包含高斯随机噪声和尖峰脉冲噪声,还含有大尺度趋势项漂移和短时波形畸变,仅靠稳健统计方法无法有效处理上述干扰.由图 12ef可见,采用经验模态分解与稳健叠加相结合的联合降噪处理,噪声干扰得到有效压制,视电阻率和视相位曲线变光滑,计算均方相对误差,分别为4.12%与2.83%, 数据质量得到明显改善.
图11 加入仿真噪声的电磁场数据(基频为32 Hz)

(a)激励电流;(b)加入仿真噪声干扰的电场;(c)加入仿真噪声干扰的磁场.

Fig 11 Electromagnetic field data with simulated noise (fundamental frequency is 32 Hz)

(a) Excitation current; (b) Electric field with simulated interferences; (c) Magnetic field with simulated noise interference.

图12 通过三种方法得到的视电阻率和视相位,并与真实值作对比

(a)和(b)为均值叠加法获得的视电阻率和相位;(c)和(d)为稳健叠加法获得的视电阻率和相位;(e)和(f)为联合降噪法获得的视电阻率和相位.

Fig 12 Apparent resistivity and apparent phase obtained by three methods and compared with the real values

(a) and (b) are the apparent resistivity and phase obtained by the mean superposition method; (c) and (d) are the apparent resistivity and phase obtained by the robust superposition method; (e) and (f) are the apparent resistivity and phase obtained by the joint noise reduction method.

综上,无论是以正演模拟生成的仿真信号作为标准信号,还是以实测高质量电磁场作为标准信号,添加噪声干扰并采用不同的算法进行降噪处理后,本文提出的联合降噪算法相比于常规的均值叠加、稳健叠加都更有优势.为了进一步测试本文信噪识别与分离算法的效果,对内蒙古某矿区约200个测点的可控源音频大地电磁数据探测进行了噪声识别和联合降噪处理.噪声识别结果表明,大部分电场数据都含有低频趋势项干扰,少数数据含有短时突发强干扰和尖峰脉冲干扰;大多数磁场数据包含尖峰脉冲干扰,少数数据包含低频趋势干扰和短时突发强干扰.通过噪声识别帮助可避免不必要的降噪处理,防止过度处理,节省时间成本.
首先选取一个测点的数据进行处理分析,测点的处理结果如图 13所示.噪声干扰主要存在于低频段.传统的均值叠加和稳健叠加都不能有效地抑制噪声干扰.联合降噪方法可以分离电磁时间序列中的噪声干扰,从而使视电阻率结果更加准确.
图13 采用三种方法对某测点实测电磁场数据进行处理所得视电阻率和视相位曲线

黑线为常规多次叠加方法,蓝线为稳健叠加方法,红线为本文联合降噪方法.

Fig 13 Apparent resistivity and apparent phase curves obtained by processing the measured electromagnetic field data of a measuring point with three methods

The black line is the conventional multiple superposition method, the blue line is the robust superposition method, and the red line is the joint noise reduction method.

最后,利用上述方法对整条测线数据进行处理,即对每个测点采集的15组基频的电场和磁场进行时域联合降噪处理,然后依次计算频域响应,再根据卡尼亚电阻率公式计算得到视电阻率和视相位的拟剖面图.为了对比不同算法的效果,同时采用均值叠加和稳健叠加两种方法对时间序列进行上述处理.图 14ab显示了均值叠加的最终计算结果.低频段(约1 Hz附近)视电阻和视相位畸变严重.图 14cd显示了通过稳健叠加方法获得的剖面.数据失真没有得到改善,反而进一步加剧.图 14ef显示了联合降噪算法的最终结果.低频段的噪声干扰被抑制,数据失真消失.总体上,经联合降噪处理后,数据质量也得到了改善,尤其是反映深部信息的低频段电阻率数据.实测数据的应用进一步验证了联合降噪算法的效果.
图14 采用三种方法对电磁场时间序列进行处理后计算得到某测线视电阻率和视相位拟断面图

(a)和(b)为常规多次叠加方法;(c)和(d)为稳健叠加方法;(e)和(f)为联合降噪方法.

Fig 14 The apparent resistivity and apparent phase pseudo profiles of a survey line are calculated after processing the electromagnetic field time series with three methods

(a) and (b) are conventional multiple stacking methods; (c) and (d) are robust stacking methods; (e) and (f) are joint noise reduction methods.

5 结论

本文提出了一种基于深度学习和时域联合降噪的综合信号处理方法,并采用仿真数据和实测数据进行了测试.首先提出了一种一维介质电磁时间序列正演算法来生成不含噪声干扰的标准参考信号.然后训练了基于长短时记忆神经网络深度学习的噪声识别器,对由低频趋势漂移、短时突发强干扰和尖峰脉冲干扰组成的八种随机混合噪声的识别精度达到95%左右,然后有针对性的调用信号降噪算法,可以避免过度处理带来的时间成本.同时,本研究还证明了基于改进经验模态分解、相关分析数据挑选和稳健统计叠加的联合降噪算法能够有效地去除上述三种典型的噪声干扰,而仅采用单一的稳健叠加方法则有可能进一步加剧信号畸变.对于仿真的含噪声干扰电磁时间序列数据,联合降噪可使数据畸变程度从20%左右降低到5%左右.对于实际探测数据,经联合降噪处理,低频段的视参数曲线更加平滑合理.

感谢中南大学、长沙巨杉科技公司等单位对本研究的帮助,感谢审稿专家对本文提出的修改意见.

Benesty J , Chen J D , Huang Y T . On the importance of the Pearson correlation coefficient in noise reduction. IEEE Transactions on Audio, Speech, and Language Processing, 2008, 16(4): 757 765

DOI

Cao J W , Cao M , Wang J Z . Urban noise recognition with convolutional neural network. Multimedia Tools and Applications, 2019, 78(20): 29021-29041

DOI

Chen C J , Jiang Q Y , Mo D . De-noising pseudo-random electromagnetic data using gray judgment criterion and rational function filtering. Chinese Journal of Geophysics, 2019, 62(10): 3854-3865

DOI

Cui Y , Liu Y , Yin C C . Denoising single-ship towed MCSEM data with adaptive frequency-division matching filter. Chinese Journal of Geophysics, 2023, 66(4): 1732-1742

DOI

Di Q Y , Wang R . Controlled Source Audio-Frequency Magneto Tellurics (in Chinese). Beijing: Science Press, 2008

Ge S C , Deng M , Chen K . Broadband signal generator for the approximation of a magnetotelluric source for indoor testing. Journal of Geophysics and Engineering, 2016, 13(4): 612-621

DOI

He J S . Wide Field Electromagnetic Method and Pseudo Random Signal Electrical Method (in Chinese). Beijing: Higher Education Press, 2010

He Z X , Chen Z C , Ren W J . Time-frequency electromagnetic (TFEM) method: Data acquisition system and its application. Oil Geophysical Prospecting, 2020, 55(5): 1131-1138

Holland P W , Welsch R E . Robust regression using iteratively reweighted least-squares. Communications in Statistics-Theory and Methods, 1977, 6(9): 813-827

DOI

Hu Y F , Li D Q , Yuan B . Application of pseudo-random frequency domain electromagnetic method in mining areas with strong interferences. Transactions of Nonferrous Metals Society of China, 2020, 30(3): 774-788

DOI

Hu Y F , Liu Z J , Li D Q . Noise separation of CSEM data based on improved clustering method. Chinese Journal of Geophysics, 2024, 67(1): 394-408

DOI

Huber P J . Robust Statistical Procedures. Philadelphia: Society for Industrial and Applied Mathematics, 1996

Jiang E H , Chen R J , Zhu D B . Static-shift suppression and anti-interference signal processing for CSAMT based on Guided Image Filtering. Earthquake Research Advances, 2022, 2(1): 100117

DOI

Jiang F B , Dong L , Dai Q W . Using wavelet packet denoising and ANFIS networks based on COSFLA optimization for electrical resistivity imaging inversion. Fuzzy Sets and Systems, 2018, 337 93-112

DOI

Key K . Is the fast Hankel transform faster than quadrature?. Geophysics, 2012, 77(3): F21-F30

DOI

Lei D , Zhao G Z , Zhang Z J . A processing method of CSAMT data with strong electromagnetc interferences by using information entropy and rational fuction filtering. Progress in Geophysics, 2010, 25(6): 2015-2023

DOI

Li G , He Z S , Deng J Z . Robust CSEM data processing by unsupervised machine learning. Journal of Applied Geophysics, 2021a, 186 104262

DOI

Li G , He Z S , Tang J T . Dictionary learning and shift-invariant sparse coding denoising for controlled-source electromagnetic data combined with complementary ensemble empirical mode decomposition. Geophysics, 2021b, 86(3): E185-E198

DOI

Li G , Liu N J , Zhang C F . Research and design of a transmission system for time-frequency-domain electromagnetic method. IEEE Access, 2019, 7 153999-154007

DOI

Li G , Wu S L , Cai H Z . IncepTCN: A new deep temporal convolutional network combined with dictionary learning for strong cultural noise elimination of controlled-source electromagnetic data. Geophysics, 2023, 88(4): E107-E122

DOI

Li Y G , Duan S M . Data preprocessing of marine controlled-source electromagnetic data. Periodical of Ocean University of China, 2014, 44(10): 106-112

Liu N . Preprocessing and research of denoising methods for marine controlled source electromagnetic data [Ph. D. thesis] (in Chinese). Changchun: Jilin University, 2015

MacLennan K , Li Y G . Denoising multicomponent CSEM data with equivalent source processing techniques. Geophysics, 2013, 78(3): E125-E135

DOI

Meier M A , Ross Z E , Ramachandran A . Reliable real-time seismic signal/noise discrimination with machine learning. Journal of Geophysical Research: Solid Earth, 2019, 124(1): 788-800

DOI

Mo D , Jiang Q Y , Li D Q . Controlled-source electromagnetic data processing based on gray system theory and robust estimation. Applied Geophysics, 2017, 14(4): 570-580

DOI

Myer D , Constable S , Key K . Broad-band waveforms and robust processing for marine CSEM surveys. Geophysical Journal International, 2011, 184(2): 689-698

DOI

Nocedal J , Öztoprak F , Waltz R A . An interior point method for nonlinear programming with infeasibility detection capabilities. Optimization Methods and Software, 2014, 29(4): 837-854

DOI

Por E , van Kooten M , Sarkovic V . Nyquist-Shannon sampling theorem. Leiden University, 2019, 1 1

Song X J , Wang X L , Dong Z . Analysis of pseudo-random sequence correlation identification parameters and anti-noise performance. Energies, 2018, 11(10): 2586

DOI

Tang J T , He J S . Controlled Source Audio Magnetotelluric Method and Application (in Chinese). Changsha: Central South University Press, 2005

Tang W W , Li Y G , Oldenburg D W . Removal of galvanic distortion effects in 3D magnetotelluric data by an equivalent source technique. Geophysics, 2018, 83(2): E95-E110

DOI

Wang M , Jin S , Wei W B . The technique analysis and achievement of the high power borehole-ground electromagnetic synchronous transmitter system. Chinese Journal of Geophysics, 2019, 62(10): 3794-3802

DOI

Wu X , Xue G Q , Zhao Y . A deep learning estimation of the earth resistivity model for the airborne transient electromagnetic observation. Journal of Geophysical Research: Solid Earth, 2022, 127(3): e2021JB023185

DOI

Xu T T , Wang Z X , Xiao Z W . Magnetotelluric power frequency interference suppression based on LSTM recurrent neural network. Progress in Geophysics, 2020, 35(5): 2016-2022

DOI

Yang Y , He J S , Li D Q . A noise evaluation method for CSEM in the frequency domain based on wavelet transform and analytic envelope. Chinese Journal of Geophysics, 2018, 61(1): 344-357

DOI

Zahn C T , Roskies R Z . Fourier descriptors for plane closed curves. IEEE Transactions on Computers, 1972, C-21(3): 269-281

DOI

Zhang B M , Jiang Q Y , Mo D . A novel method for handling gross errors in electromagnetic prospecting data. Chinese Journal of Geophysics, 2015, 58(6): 2087-2102

DOI

Zhang X , Li D Q , Li J . Signal-noise identification for wide field electromagnetic method data using multi-domain features and IGWO-SVM. Fractal and Fractional, 2022, 6(2): 80

DOI

Zhao H T , Liu L H , Wu K . Constant voltage-clamping bipolar pulse current source for transient electromagnetic system. Electric Power Components and Systems, 2013, 41(10): 960-971

DOI

Zhou C , Tang J T , Pang C . A theory and simulation study on the space-time array hybrid source electromagnetic method. Chinese Journal of Geophysics, 2019, 62(10): 3827-3842

DOI

Zhou H G , Yao Y , Liu C S . Feasibility of signal enhancement with multiple grounded-wire sources for a frequency-domain electromagnetic survey. Geophysical Prospecting, 2018, 66(4): 818-832

DOI

Ziolkowski A , Wright D , Mattsson J . Comparison of pseudo-random binary sequence and square-wave transient controlled-source electromagnetic data over the Peon gas discovery, Norway. Geophysical Prospecting, 2011, 59(6): 1114-1131

DOI

超健 , 奇云 , . 基于灰色判别准则和有理函数滤波的伪随机电磁数据去噪. 地球物理学报, 2019, 62(10): 3854-3865

DOI

, , 长春 . 单船拖曳可控源电磁数据自适应分频匹配滤波消噪方法. 地球物理学报, 2023, 66(4): 1732-1742

DOI

青云 , . 可控源音频大地电磁数据正反演及方法应用. 北京: 科学出版社, 2008

继善 . 广域电磁法和伪随机信号电法. 北京: 高等教育出版社, 2010

展翔 , 忠昌 , 文静 . 时频电磁(TFEM)勘探技术: 数据采集系统. 石油地球物理勘探, 2020, 55(5): 1131-1138

艳芳 , 子杰 , 帝铨 . 基于优化聚类的人工源电磁法数据信噪分离方法. 地球物理学报, 2024, 67(1): 394-408

DOI

, 国泽 , 忠杰 . 强干扰地区CSAMT数据信息熵与有理函数滤波的处理方法. 地球物理学进展, 2010, 25(6): 2015-2023

DOI

予国 , 双敏 . 海洋可控源电磁数据预处理方法研究. 中国海洋大学学报(自然科学版), 2014, 44(10): 106-112

. 海洋可控源电磁数据典型预处理及几种去噪方法研究[博士论文]. 长春: 吉林大学, 2015

井田 , 继善 . 可控源音频大地电磁法及其应用. 长沙: 中南大学出版社, 2005

, , 文博 . 大功率井-地电磁同步发射技术分析与系统实现. 地球物理学报, 2019, 62(10): 3794-3802

DOI

滔滔 , 中兴 , 卓伟 . 基于LSTM循环神经网络的大地电磁工频干扰压制. 地球物理学进展, 2020, 35(5): 2016-2022

DOI

, 继善 , 帝铨 . 在频率域基于小波变换和Hilbert解析包络的CSEM噪声评价. 地球物理学报, 2018, 61(1): 344-357

DOI

必明 , 奇云 , . 电磁勘探数据粗大误差处理的一种新方法. 地球物理学报, 2015, 58(6): 2087-2102

DOI

, 井田 , . 时空阵列混场源电磁法理论及模拟研究. 地球物理学报, 2019, 62(10): 3827-3842

DOI

Outlines

/