组合广义形态滤波在鹿井地区音频大地电磁测量数据去噪中的应用

2022-01-27李于波刘祜段书新汪硕

铀矿地质 2022年1期

李于波，刘祜，段书新，汪硕

（核工业北京地质研究院中核集团铀资源勘查与评价技术重点实验室，北京 100029）

鹿井地区位于湖南省郴州市汝城县境内，处于湖南、江西、广东三省交界处，具有铁、钨、铜、稀土、铀等多种矿产资源。为查明鹿井铀矿田深部成矿要素，开展了音频大地电磁（Audio Magnetotelluric Method，简称AMT）测量。但由于工作区内电磁干扰严重，部分采集到的时间序列上叠加了与地质体无关的噪音信号，需进行去噪处理以便反演处理得到地下的真实电性结构。

针对音频大地电磁去噪，前人从时间域或频率域的角度先后开展了互功率谱法、小波分析法、远参考技术、小波自适应阈值、改进阈值的TI 小波、时间序列依赖关系等去噪技术［1-7］。上述系列去噪方法取得了一定的应用效果，但亦存在若干问题，如：互功率谱法只能去除不相关噪声、小波变换参数选取不当会导致去噪效果适得其反等。为实现对大量实测数据的快速去噪，笔者以鹿井地区多个测点的时间序列为研究对象，探讨了组合广义形态滤波方法在该区AMT 去噪中的应用。

1 鹿井地区AMT 工程部署及电磁干扰特征

1.1 AMT 工程部署

为查明鹿井铀矿田深部地质情况，在该区布设了4 条AMT 探测剖面（图1）。野外利用V8（主机）结合AMTC30（磁探头）的方式开展数据采集，获取到TS2（高频）、TS3（中频）、TS4（低频）3 个时间序列文件，对应的采样率分别为24 000、2 400、150 Hz。

图1 鹿井地区AMT 测线布置图Fig.1 Layout of AMT survey line in Lujing area

1.2 AMT 电磁干扰特征

据现场实际踏勘，测区内有高速、省道、乡道等公路贯穿，同时分布了较多的广播电台、高压线、发射塔等，预计会对AMT 实测信号产生较多干扰。对实测的时间序列进行分析，证实了上述干扰的存在。

图2a 中，Ex电场信号由于受到工频噪声干扰，其背景值明显高于其他测点的背景水平。图2b 中，Ey电场信号在180 和400 采样点处出现明显的尖波，为三角波噪音在时间序列上的反映。图2c 中，Hx磁场信号在30 和330 采样点处亦出现2 个格值超过150 000 的脉冲。图2d中，Hx磁场信号出现规律性的似充放电三角波噪声干扰。为获取地下真实的电性结构，对上述含噪音数据进行去噪处理就显得尤为重要。

图2 鹿井地区采集电场信号受多种噪声干扰图Fig.2 The acquired signal interfered by various noises in Lujing area

2 组合广义形态滤波原理及其在AMT去噪中的具体实现

形态滤波是由法国数学家Matheron G和Serra J于1964年共同创立的信号分析方法，其核心思想是建立适当的结构元素，让待处理的序列（或集合）依次穿过结构元素，从而提取到该序列（或集合）的数学形态。在数据分析过程中，若同时使用不同的结构元素进行组合处理，即为组合广义形态滤波。

在AMT 时间序列处理过程中，为克服时间序列的基线漂移，采用正负结构元素串行的方法构建组合广义形态滤波器［8］，其流程如图3 所示。以Ex 电场信号为例，将信号Input_Ex 依次通过ΨGOC(GCO)(s1，s2)和ΨGOC(GCO)(-s1，-s2)两个基本滤波单元（其中S1，S2 为结构元素），即得到Ex 电场中所包含的噪声序列Exlb(n)。

图3 组合广义形态滤波流程图Fig.3 Flow chart of combined generalized morphology filtering

图（3）中，ΨGOC(GCO)(s1，s2) 滤波单元的定义为

其中GOC(f(n))、GCO(f(n))分别为广义形态开-闭滤波器和广义形态闭-开滤波器，“∘”为形态开运算，“•”为形态闭运算，其定义如式（2）、（3）所示：

式（2）和式（3）中的Θ 为腐蚀运算，其目的是减少待处理序列的峰值，即剔除掉不平滑的峰值数据，缩小峰域和拓宽谷域，使目标序列的值收缩；⊕为膨胀运算，其目的是拓宽峰域和缩小谷域，使目标序列的值膨胀，两运算的公式分别如式（4）、（5）所示：

在本次AMT 去噪处理中，将S1 和S2 分别设定为5 点圆盘型结构元素（图4a）和5 点抛物线型结构元素（图4b），两者的函数表达式分别如式（6）、（7）所示：

图4 五点圆盘型（a）和五点抛物线型（b）结构元素示意图Fig.4 Diagram of structuring elements of disc type（a）and parabolic type（b）

据图3 获取到噪声序列Exlb(n)后，将其从Ex电场信号中剔除，即得到重构后的时间序列信号Exdn(n)。

3 实测数据的去噪应用

为实现对大量实测数据的快速去噪，笔者以Python 编程语言为平台，按照图5 所示流程实现了鹿井地区AMT 实测数据的组合广义形态滤波去噪。整个流程的核心是对4个电磁场分量做正、负结构元素的广义形态滤波处理，并重构TS 文件。其中TS2 时间序列的结构元素参数L2=3，k2=0.5；TS3 的结构元素参数L3=5，k3=0.4；TS4 的结构元素参数L4=2，k4=0.3。

图5 数据处理及组合广义形态滤波去噪流程图Fig.5 Flow chart of data processing and filtering denoising by generalized morphology combination

以line2-05、line2-07、line2-27 测点为例，重构后的时间序列（图6、7、8）分别去除了阶跃噪声、三角波、似充放电三角波的干扰，表明组合广义形态滤波对该类电磁干扰具有较好的去噪效果。

图6 鹿井地区采集的含噪时间序列去噪前、去噪后信号对比图（line2-05）Fig.6 The comparison of original series and denoised series in Lujing area（line2-05）

图7 鹿井地区采集的含噪时间序列去噪前、去噪后信号对比图（line2-07）Fig.7 The comparison of original series and denoised series in Lujing area（line2-07）

图8 鹿井地区采集的含噪时间序列去噪前、去噪后信号对比图（line2-27）Fig.8 The comparison of original series and denoised series in Lujing area（line2-27）

上述去噪效果在各个测点的功率谱密度上亦有所表现。以Line2-07 测点为例，图9a 展示了TS2 序列中各分量的功率谱密度，曲线在50 Hz 及其倍频处出现了明显的峰值情况，表明该点电磁场信号受到了工业谐波干扰；同时，低频的功率谱密度小、高频的功率谱密度大，推测是三角波等大格值噪声淹没了低频有用信号所致。对其进行组合广义形态滤波处理，得到重构后的时间序列功率谱密度曲线如图9b 所示。可以看出，去噪后的时间序列在150、250、350、450 Hz 等频率处的工业噪声得到了压制，且低频1～300 Hz 范围内的功率谱密度得到了加强，Ex、Ey、HxHy 四道信号的功率谱密度均较为稳定，分别保持在1 000、100、10 000、1 000 w/Hz 左右。

图9 Line2-07 测点TS2、TS3 序列去噪前、后功率谱密度对比图Fig.9 The comparison of original and denoised TS2 and TS3 PSD of Line 2-07 point

相比较而言，该测点处TS3 时间序列受噪声影响较小，经组合广义形态滤波处理前后的功率谱密度曲线变化不大。

图10 为测点Line2-05 去噪前后的视电阻率曲线对比。由图可知：在小于1 Hz 的频率范围，rhoxy 和rhoyx 经广义形态滤波去噪后变得更加光滑、稳定。针对去噪后rhoxy 曲线在10 Hz 左右出现的跳变，推测是由于形态滤波时将部分有用信号进行了剔除所致，需要在后续数据处理过程中修改滤波参数，以减少有效信号的损失，进一步提高广义形态滤波效果。