基于自适应多分辨率奇异值分解的大地电磁数据处理

李晋，马翻红，汤井田，李勇

1 湖南师范大学信息科学与工程学院，长沙 410081 2 中南大学有色金属成矿预测与地质环境监测教育部重点实验室，长沙 410083 3 中国地质科学院地球物理地球化学勘查研究所，自然资源部地球物理电磁法探测技术重点实验室，河北廊坊 065000

0 引言

20世纪50年代，Tikhonov和Cagniard提出了大地电磁测深法（Magnetotelluric,MT）（Tikhonov,1950;Cagniard,1953）；该方法是一种探测地壳深部结构的地球物理方法（赵国泽等,2004），具有探测深度大、成本低、施工方便、垂向分辨能力和水平分辨能力高等优点，在地震预测、油气田勘探与普查、矿产资源勘查等领域已得到广泛应用（汤吉等,2005;魏文博等,2009;金胜等,2010;胡祥云等,2012;董树文等,2012）.然而，频带范围极宽的天然大地电磁信号非常微弱，且不可避免地会受到各类噪声的干扰，尤其是在矿集区，大地电磁数据质量急剧下降，导致视电阻率-相位曲线发生畸变、阻抗估计过度失真（白登海等,2002）.此时的数据已不能客观反映地下电性结构，严重影响了后续电磁反演成像的可靠性和可解释性，给地球物理勘探带来极大困扰（汤井田等,2012b;吕庆田等,2015）.因此，有效压制强干扰，提升数据质量是大地电磁测深领域最具挑战的任务之一，对精细勘探地下深部结构、勘查深部矿产资源具有重要的理论意义和实际应用价值.

国内外学者针对不同的噪声提出了许多有针对性的去噪方法.Gamble等（1979）提出了远参考大地电磁测深法，该方法用于去除同源相关电磁噪声；Egbert和Booker（1986）提出了消除非高斯分布噪声的Robust统计方法；Kao和Rankin（1977）提出了利用最小二乘法中的循环计算来消除非相关噪声；由于小波变换的多分辨特性，Trad和Travassos（2000）将其用于非平稳的大地电磁数据处理，随后该方法在阻抗估计和长（短）周期噪声压制等方面均得到广泛应用（徐义贤和王家映,2000;范翠松等,2008;凌振宝等,2016;Ling et al.,2019）；王书明和王家映（2004）将高阶统计量运用于大地电磁数据处理，有效抑制了高斯有色噪声；为了压制工频干扰，汤井田等（2008）和Cai（2017）将Hilbert-Huang变换引入大地电磁数据处理；紧接着，互补集成经验模态分解等系列方法在大地电磁数据处理中得以应用，该方法解决了集成经验模态分解（史恒等,2011）出现的分解误差，以及经验模态分解（Empirical Mode Decomposition,EMD）出现的模态混叠和端点效应（Li et al.,2019）；汤井田等利用数学形态滤波有效抑制了具有明显形态特征的大尺度干扰（汤井田等,2012a），并结合稀疏表示和字典学习对大地电磁数据进行去噪（汤井田等,2018）；李晋等研究了匹配追踪和变分模态分解（Variational Mode Decomposition,VMD）的大地电磁信噪分离方法（李晋等,2018,2019;Zhang et al.,2021）.

分析国内外相关文献可知，大地电磁噪声强度大且环境复杂，尤其是在强干扰区，现有方法具有一定的优势，同时也存在一定的局限性.

近年来，奇异值分解（Singular Value Decomposition,SVD）在信号处理领域迅猛发展.与小波分析相比，该方法具有稳健性强、信噪比高、零相移、波形失真小等优点，在故障诊断、电能质量检测、信号识别和消噪，以及数据挖掘等诸多领域得到广泛应用.Phillips等（2009）将SVD应用于卫星遥感数据降维，获得了比主成分分析（王文波等,2013）更精确的特征，分类速度更快且更准确；Lehtola等（2008）将数学形态学和SVD相结合，对心电信号进行消噪处理，提高了心电信号的准确性.赵学智等（2010,2017）在故障诊断领域，结合SVD对矩阵分解及相空间矩阵构造等进行探索，并将小波的多分辨思想应用于SVD，提出了多分辨率奇异值分解（Multi-Resolution Singular Value Decomposition,MRSVD），同时分析了SVD和小波变换在信号处理方面的优势.Luo和Zhang（2019）将MRSVD应用于滚动轴承周期脉冲特征提取，实现了初期的故障检测.He等（2020）采用MRSVD抑制振动噪声，准确地检测了振动信号中的故障点，并利用长短时记忆神经网络预测了轴承性能；Zhang等（2020）研究了各种多分结构的奇异值去噪模型.鉴于MRSVD在去噪方面的优越性，本文结合MRSVD和标准差差值提出基于自适应多分辨率奇异值分解（Adaptive Multi-Resolution Singular Value Decomposition,AMRSVD）的大地电磁数据处理方法，并与EMD、VMD等方法进行分析对比.结果表明，该方法不受噪声先验信息的影响、处理效率高，且能有效分离相关性较好的噪声，保留了大地电磁数据的原始特征.

1 算法原理

1.1 奇异值分解

Kaleman（1996）提出奇异值分解，该方法原理描述如下：

对任意一个矩阵H∈Rm×n，都有正交矩阵U=（u1,u2,…,um）∈Rm×m及正交矩阵V=（v1,v2,…,vn）∈Rn×n，使得矩阵H满足

H=UΣVT，

（1）

式中，Σ=（diag（α1,α2,…,αr）,0）（Σ∈Rm×n），r=min（m,n），其中α1≥α2≥…≥αr>0，αi（i=1,2,…,r）为矩阵H的奇异值.

对信号进行SVD处理时，首先将信号构造Hankel矩阵H.设待处理信号x（t）=[x1,x2,x3,…,xN]，矩阵H的形式如下：

（2）

式中1

一般地，只对信号构造一次Hankel矩阵，然后运用SVD进行处理.从本质上讲，这种对信号采用单一的SVD处理，其结果是利用正交矩阵U（V）的行（列）矢量得到的，都属于同层次的矢量空间.

1.2 多分辨率奇异值分解

为了将SVD拓展到不同层次的矢量空间，赵学智等（2010）借鉴小波的多分辨率特征，在SVD的基础上提出一种二分递推矩阵构造方法;Zhang等（2020）在二分递推矩阵的基础上提出多分递推矩阵构造方法，用于将复杂信号分解到不同层次的子空间.

1.2.1 三分递推SVD算法

三分递推SVD算法步骤如下：

（1）对待处理数据x（t）用三分递推构造一个Hankel矩阵H1：

（3）

（2）对矩阵H1进行SVD处理：

（4）

式中，U1=|_u1,1,u1,2,u1,3_|（U1∈R3×3）、V1=|_v1,1,v1,2,v1,3_|（V1∈R（N-2）×3）分别为对应的左正交矩阵和右正交矩阵，相应的对角矩阵为Σ1=（diag（α1,1,α1,2,α1,3）,0）（Σ1∈R3×3）.其中，α1,1、α1,2、α1,3分别为信号第1次分解得到的奇异值.

（3）通过式（4）得到三个奇异值α1,1、α1,2、α1,3，并进行SVD重构得到：

（5）

式中，u1,1、u1,2和u1,3为3×1的矩阵，分别表示u1,1=（u1,1,1,u1,1,2,u1,1,3）T、u1,2=（u1,2,1,u1,2,2,u1,2,3）T及u1,3=（u1,3,1,u1,3,2,u1,3,3）T，将其代入式（5）得到：

（6）

同理，得到HD1和Hd1.HA1对应的奇异值大，反映信号的主体成分，为近似矩阵；HD1和Hd1对应的奇异值小，反映信号的细节成分，为细节矩阵；从这三个矩阵可以分别恢复出近似信号A1和细节信号D1、d1.

（4）通过上述步骤得到第1次分解的结果，再利用近似信号A1构造Hankel矩阵H2返回步骤（2）进行同样的处理，可将原始信号分解为不同分辨率的近似信号和细节信号，如图1所示.

图1 三分递推SVD分解流程Fig.1 The decomposition process of three recursive SVD

1.2.2 三分递推SVD的多分辨率特性

多分辨率指的是L2（R）的子空间Vj中每个子空间Vj+1都是它前一级Vj的精细化.从式（6）可知，在三分递推SVD中，近似信号Aj是利用vj,1得到的，而细节信号Dj和dj分别是由vj,2和vj,3得到；vj,1、vj,2及vj,3称为近似基矢量和细节基矢量，通过这些基矢量可得到不同分辨率的近似信号和细节信号.从三分递推SVD的分解过程可知，下一层分解的基矢量是由上一层的近似基矢量得到，说明三分递推SVD在分解过程中实现了类似于小波的多分辨率分解，即多分辨率奇异值分解（MRSVD）.

1.2.3 三分递推SVD的分解性能

为了测试MRSVD的优越性，构造如图2所示的原始信号做模拟仿真实验，其中f1=sin（80πt），f2=sin（300πt）×（1+2sin（30πt）），f3=f1+f2；左边为时域波形图，右边为其对应的频谱.

将f3（合成信号）作为分解对象，采用EMD、VMD、SVD和MRSVD四种方法进行处理，如图3所示；其中左边为时域波形图，右边为其对应的频谱.分析图3可知，（a）中EMD分解得到的IMF1和IMF2出现端点效应；（b）中VMD虽克服了EMD分解存在的缺陷，但单个分量的频率成分不连续；（c）中SVD将原始信号分成了三个分量，每个分量虽都包含原信号信息，但仍无法得出精确的分解效果；（d）中MRSVD将原始信号分解为A1～A7七个近似信号，最后一个分量为细节信号的叠加；原始信号经MRSVD分解得到的近似信号A6的时间域波形与f1的相似度达到0.9921，A7与f1的相似度为0.9927，且（d）中最后一个分量（f3-A7）与f2的时间域波形相似度达到0.9866，说明信号分解到A7时对应的细节信号的频谱即f3-A7与f2最为接近.分析不同方法的分解效果可知，MRSVD能将原始信号中的复杂信号以细节信号的形式分解出来，凸显了隐藏的原始信号特征.

1.3 标准差

标准差在概率统计中通常用来作为统计分布程度上的测量，反映数据集个体间的离散程度.标准差越大，表示数据和平均值之间的差异较大；反之，则表示越接近平均值.

（7）

2 模拟仿真实验

2.1 大地电磁信号分布特性

大地电磁场指的是由地球外部场源引起天然电磁场短周期变化形成的电场和磁场，其频率范围极

图2 原始信号的时域波形和频谱Fig.2 Time domain waveform and frequency spectrum of original signal

（续图3）

图3 不同分解方法的时域波形和频谱对比图（a） EMD;（b） VMD;（c） SVD;（d） MRSVD.Fig.3 Comparison of time domain waveform and frequency spectrum of different decomposition methods

宽、信号极其微弱，且场源极化方向随机；同时具有复杂的时变特性，野外采集的大地电磁信号极易受到各种电磁噪声的污染.由于天然大地电磁信号是典型的非线性、非平稳信号，实测数据受到的大尺度噪声和低频有用信号在统计分布上具有明显区别.

图4所示为一组测试样本，分别为青海地区某测点中一段几乎未受电磁干扰的大地电磁数据段和矿集区测点中含类充放电三角波干扰、类脉冲干扰、类方波干扰的数据段.

图5所示为图4中的测试样本经式（3）和（4）得到的三个奇异值α1,1、α1,2、α1,3的分布特性，其中每200个采样点为1段，共分为15段.分析图5可知，无干扰信号的三个奇异值在数值上没有明显差异；类充放电三角波干扰α1,1明显大于α1,2和α1,3，且α1,2和α1,3非常接近；类方波干扰和类脉冲干扰在有干扰的部分，α1,1也明显大于α1,2和α1,3；同样地，图6所示为青海点QH401504中一段电道Ex和磁道Hx数据经式（3）和（4）得到的三个奇异值α1,1、α1,2、α1,3的分布特性，其中左边为时域波形图，右边为其对应的奇异值分布特性.从图6可知，电道和磁道数据在受干扰的情况下同图5相比，其奇异值分布具有相似性.分析图5和图6可知，奇异值大的可表征干扰数据，奇异值小的则可代表低频信号成分.由于α1,2和α1,3分别为Dj和dj的奇异值，且相差不大，为此后续将细节信号定义为Dj和dj的叠加.

图7所示为图4测试样本中的无干扰信号、类充放电三角波干扰、类脉冲干扰和类方波干扰数据分别采用MRSVD分解1层得到的近似信号和细节信号.

分析图7可知，图4中的无干扰信号经MRSVD得到的细节信号和近似信号没有明显区别；类充放电三角波干扰分解得到的近似信号主要表现为强干扰，细节信号则主要由微弱的有用信号构成；类脉冲和类方波干扰分解之后的细节信号中出现大量的脉冲干扰，分解效果不佳.

为了区分强干扰和微弱的大地电磁有用信号，分别计算图7中的近似信号标准差与细节信号标准差的差值（简称标准差差值），如图8所示.

分析图8可知，经MRSVD分解1层后，无干扰信号的标准差差值在基线附近，类充放电三角波干扰基本都在0.5以上；类脉冲干扰在采样点为500时有明显干扰，但差值仍在基线附近，干扰数据区分不明显；类方波干扰在采样点0～500处有明显干扰，其标准差差值也在0.5以上.为此，文中采用MRSVD分解1层得到的标准差差值对类充放电三角波干扰和类谐波干扰进行信噪辨识，其辨识参数定义为：

图4 测试样本数据Fig.4 Test sample data

图5 测试样本奇异值分布特性（a） 无干扰；（b） 类充放电三角波；（c） 类脉冲；（d） 类方波.Fig.5 Singular value distribution characteristics of test sample data（a） No interference;（b） Charge-discharge-like triangle wave;（c） Pulse-like wave;（d） Square-like wave.

图6 测点QH401504中数据段的奇异值分布特性Fig.6 Singular value distribution characteristics of data segments in site QH401504

图7 经MRSVD分解1层得到的近似信号和细节信号Fig.7 The approximate signal and detail signal obtained by MRSVD decomposition layer 1

|Δψ|=|ψA1-ψ（D1+d1）|，

（8）

约束条件定义为|Δψ|<θ（θ为0.4～0.9），其中ψA1、ψ（D1+d1）分别为近似信号标准差和细节信号标准差.

图9所示为两段实测数据（含类充放电三角波干扰和类谐波干扰）采用上述方法和约束条件得到的信噪辨识效果.分析图9可知，有用信号段和强干扰数据段可以有效区分.

图8 近似信号标准差和细节信号标准差的差值对比Fig.8 Comparison of difference between approximate signal standard deviation and detail signal standard deviation

图10所示为测试样本（图4）中含类充放电三角波干扰的原始数据经MRSVD分解11层得到的近似信号.分析图10可知，类充放电三角波干扰的噪声轮廓随着分解层数的增大逐渐分解到近似信号部分.分析A1～A8可知，噪声轮廓逐渐光滑，低频信号被分解到细节信号部分，A8～A11的噪声轮廓则没有发生明显变化，说明一定的分解层数对MRSVD的去噪效果敏感.

图11所示为图9中两段实测数据经MRSVD分解50层得到的相邻细节信号标准差差值的变化情况.

分析图11可知，两段含强干扰数据的相邻细节信号标准差差值随分解层数的增加，差值增大随后逐渐减小直到0.001～0.01趋于稳定状态.结合图10分析可知，噪声轮廓随着分解层数的增加逐步变得光滑，但当分解到一定层数后，近似信号不再变化，对应的细节信号也不再变化.考虑到大地电磁信号的主要特征集中在细节信号（微弱的低频成分），为此文中将相邻细节信号标准差差值|Δψ（Dj+dj）|作为最佳分解层数的约束条件（刘嫣和汤伟,2016），即若存在分解层数j使得：

图9 不同干扰的辨识效果Fig.9 Identification effect of different interference

|Δψ（Dj+dj）|<ω，

（9）

则MRSVD停止分解，即完成自适应多分辨率奇异值分解（AMRSVD），式中ω为0.001～0.01.

2.2 算法流程

基于AMRSVD的大地电磁数据处理算法流程如图12所示.

具体步骤如下：

（1）初始化参数：分解层数j=1、ω取0.001～0.01、θ取0.4～0.9，对大地电磁数据均匀分段；

（2）对大地电磁数据构建Hankel矩阵H，对其进行SVD分解得到近似信号和细节信号；

（3）计算近似信号和细节信号的标准差差值|Δψ|，判断是否满足|Δψ|<θ，若满足则认为是大地电磁有用信号，否则进入步骤（4）进行信噪分离处理；

图10 原始信号经MRSVD分解得到的近似信号Fig.10 The approximate signal obtained by MRSVD decomposition of the original signal

图11 相邻细节信号标准差差值变化Fig.11 Variation of standard deviations of adjacent detail signal

（4）对分解出的近似信号重建Hankel矩阵且j=j+1，利用SVD分解得到下一层的近似信号和细节信号；

（5）计算相邻细节信号的标准差差值|Δψ（Dj+dj）|，判断是否满足|Δψ（Dj+dj）|<ω，若满足则停止分解，得到最终的近似信号即噪声轮廓，然后用原始数据减去近似信号得到低频信号，并与步骤（3）中的有用信号进行拼接、重构，否则返回步骤（4）.

2.3 EMTF数据模拟

将EMTF开源代码包中的纯净电道（Ex、Ey）数据和磁道（Hx、Hy）数据分别加入相关性较强的大尺度噪声进行分析，引入均方误差（MSE）、信噪比（SNR）、归一化互相关性（NCC）和运行效率（Runtime）进行评判（Li et al,2020a,2020b,2021）.

图13和图14所示分别为EMTF含噪数据（Ex、Ey）采用EMD、VMD、SVD和AMRSVD的时频域去噪效果，其中左边为时域波形图，右边为其对应的频谱.分析可知，经EMD处理后，时域波形中大部分有用信号丢失严重、频谱失真；经VMD和SVD处理后残留了大量的尖脉冲，去噪性能较低；AMRSVD可以较好地实现大尺度强噪声和低频信号的分离.

图15所示为采用不同方法处理图13和图14的含噪EMTF数据得到的视电阻率曲线对比图.分析图15可知，原始EMTF数据的视电阻率曲线呈平滑直线状态，当Ex和Ey分别加入大尺度噪声后视电阻率曲线发生了显著跳变.经EMD和SVD处理后，视电阻率曲线呈下降趋势，低频有用信号丢失；VMD虽然比EMD和SVD有较大提升，但中低频段处理不佳；相比而言，AMRSVD得到的视电阻率曲线改善明显.

图12 AMRSVD数据处理流程图Fig.12 Flow chart of AMRSVD data processing

图13 Ex道数据经不同方法处理的时频域去噪效果Fig.13 Denoising effect of different methods in time-frequency domain of Ex data

图15 EMTF数据的视电阻率曲线对比图Fig.15 Comparison of apparent resistivity curve for EMTF data

图16 实测数据时间域去噪效果图（a） 类充放电三角波；（b） 类谐波.Fig.16 Time domain denoising effect of measured data（a） Charge-discharge-like triangle wave;（b） Harmonic-like wave.

图17 四道实测数据段去噪前后对比（a） 去噪前；（b） 去噪后.Fig.17 Comparison of four measured data segments before and after denoising（a） Before denoising；（b） After denoising.

图18 测点BL22200J的视电阻率-相位曲线对比图Fig.18 Comparison of apparent resistivity-phase curve for site BL22200J

图19 测点EL22189A的视电阻率-相位曲线对比图Fig.19 Comparison of apparent resistivity-phase curve for site EL22189A

图20 测点EL22195A的视电阻率-相位曲线对比图Fig.20 Comparison of apparent resistivity-phase curve for site EL22195A

表1和表2分别为Ex道和Ey道数据经不同方法的去噪性能对比.分析表1、表2可知，AMRSVD在NCC、MSE、SNR均优于EMD、VMD和SVD的去噪效果，尤其是AMRSVD的去噪速度更快、效率更高.

表1 Ex道数据不同方法的去噪性能对比Table 1 Comparison of denoising performance of different methods of Ex data

表2 Ey道数据不同方法的去噪性能对比Table 2 Comparison of denoising performance of different methods of Ey data

3 实测数据处理

3.1 时间序列分析

为了验证方法的实用性，将其应用于庐枞矿集区相关性较强的大尺度类充放电三角波干扰和类谐波干扰的大地电磁时间序列，如图16所示.

分析图16a可知，AMRSVD可以有效去除大地电磁数据中的类充放电三角波干扰、保留更多的有用信号.分析图16b可知，AMRSVD对一些类谐波干扰也可以提取出光滑的噪声轮廓曲线.

图17所示为庐枞矿集区某测点中同一时段的电道和磁道数据经AMRSVD法去噪前后的时域波形对比图.

分析图17a可知，Ex/Hy和Ey/Hx表现出很强的噪声相关性；分析图17b可知，四道数据中的强干扰均被有效去除，低频信号得到了更好保留.

3.2 视电阻率-相位曲线分析

图18—20所示分别为庐枞矿集区测点BL22200J、EL22189A和EL22195A经远参考（RR）、EMD、VMD、SVD和AMRSVD处理的视电阻率-相位曲线对比图.

分析图18—20可知，3个测点的原始数据视电阻率曲线在低频段急剧上升，30～0.5 Hz的视电阻率曲线基本呈45°，视电阻率值从102Ωm上升到了106Ωm，提升了4个数量级，对应的相位在0°或-180°附近且分布紊乱，这些测点的数据真实性急剧下降，表现为典型的近源效应.经RR处理后，视电阻率曲线在30～5 Hz变得连续、光滑，但在5～0.3 Hz处的视电阻率曲线跳跃大且分布凌乱，表明远参考对近源干扰无能为力.经EMD处理后，视电阻率曲线在30～5 Hz较为连续，但小于5 Hz的曲线部分频点的值急剧下降且相位较为紊乱，低频信息丢失严重.经VMD处理后，30～5 Hz视电阻曲线仍呈45°上升，3～0.3 Hz视电阻率曲线呈下降趋势，去噪效果不佳.SVD由于分解矩阵难以确定，导致SVD的去噪效果较差，视电阻率曲线紊乱.分析庐枞矿集区大量原始大地电磁数据可知，电道中含有大量的类充放电三角波噪声，磁道中含有不同程度的类充放电三角波噪声和类谐波噪声，且相关性强，这些噪声可能是引起近源效应的主要原因，经AMRSVD处理后视电阻率曲线45°上升的趋势得到明显缓解.然而由于该矿集区采集的原始大地电磁数据中广泛分布脉冲干扰，导致1～10 Hz左右的频段数据严重受损，而文中方法对脉冲噪声分离效果不佳，重构的信号中残留有尖脉冲干扰.

4 结论

针对如何从强干扰中有效提取微弱的大地电磁信号，文中以类充放电三角波干扰和类谐波干扰作为研究对象，提出一种基于自适应多分辨率奇异值分解的大地电磁数据处理方法.研究了大地电磁数据的奇异值分布特性，发现奇异值的大小能分别表征近似信号和细节信号；通过将大地电磁数据分解成不同分辨率的近似信号和细节信号，并引入标准差差值对大地电磁数据进行信噪辨识；同时，结合相邻标准差差值和MRSVD，提出利用AMRSVD对强干扰数据进行去噪处理.仿真实验和实测数据处理结果表明，本文方法能有效去除时间序列中形状规则、相关性较强的干扰，视电阻率-相位曲线得到改善，去噪后的信号可靠性提升；与EMD、VMD和SVD等方法相比，AMRSVD的处理效果更优、分解效率更高，为强干扰下开展大地电磁数据处理提供了一种新的研究思路.

由于本文方法的处理效果与噪声的相关性及出现的频率等有关，若噪声的能量、幅值等特征不明显，辨识度将降低、信噪分离会有一定的偏差.同时，该方法对尖脉冲干扰处理效果不佳.为此，接下来将重点研究方法对不同噪声类型的稳健性，以及引入智能算法优化MRSVD的自适应分解过程，提升方法的抗噪性能.