高频大地电磁噪声压制

引用本文

化希瑞, 席振铢, 胡双贵, 肖晓. 高频大地电磁噪声压制[J]. 物探与化探, 2017,41(2): 328-332.
HUA Xi-Rui, XI Zhen-Zhu, HU Shuang-Gui, XIAO Xiao. Noise suppression of high frequency magnetotelluric method[J]. Geophysical and Geochemical Exploration, 2017,41(2): 328-332. 复制到剪切板

Doi:10.11720/wtyht.2017.2.21
Permissions

《物探与化探》编辑部所有

高频大地电磁噪声压制

化希瑞^1,², 席振铢¹, 胡双贵¹, 肖晓¹

1.中南大学地球科学与信息物理学院, 湖南长沙 410083

2.中铁第四勘察设计院集团有限公司,湖北武汉 430063

作者简介: 化希瑞(1982-),博士研究生,研究方向为电磁法信号处理与数据处理。Email:huaxirui@qq.com

收稿日期: 2016-12-12

基金: 国家高技术研究发展计划(“863”计划)项目(2014AA06A602); 中国铁建股份有限公司科技研究开发计划课题(14-C61)

摘要

高频大地电磁法(HMT)信号在20 kHz~100 kHz频率范围内易受人文噪声干扰,电磁信号强度弱,数据误差大。本文提出采用高阶谱分析重构高频大地电磁非最小相位信号,进而有效抑制随机高斯噪声、提高信噪比。通过数值模拟信号和实测数据处理,证明该方法能有效压制噪声,在提取实际地质信息方面优于传统方法。

关键词: 高频大地电磁法; 高阶谱; 信号处理; 噪声压制

中图分类号:P631 文献标志码:A 文章编号:1000-8918(2017)02-0328-05

Noise suppression of high frequency magnetotelluric method

HUA Xi-Rui^1,², XI Zhen-Zhu¹, HU Shuang-Gui¹, XIAO Xiao¹

1.School of Geosciences and info-Physics, Central South University, Changsha 410083, China

2. China Railway Siyuan Survey and Design Group Co., Ltd., Wuhan 430063, China

Abstract

In the frequency range of 20kHz ~ 100kHz, the HMT signal has the characteristics of weak strength, low correlation and poor repeatability, which is caused by the interference of human noise in the near surface. In this paper, the authors used higher-order spectrum method to analyze and reconstruct HMT signal, which can effectively suppress Gaussian random noise and raise signal noise ratio. Numerical analog signal and measured data processing prove that the method is feasible. According to the application effect, the method is better than the traditional method in colored noise suppression and practical geological information extraction.

Keyword: high frequency magnetotelluric method; higher-order spectrum; signal processing; noise suppression

Show Figures

0 引言

高频大地电磁法(HMT)是根据深埋隧道的工程地质勘探需要提出来的, 工程勘察最关心的是近地表1 000 m内的地层地质情况^{[1, 2]}, 要求大地电磁法观测频段为10 Hz~100 kHz, 与传统的大地电磁法(MT)观测频段为0.001~340 Hz、音频大地电磁法(AMT)的观测频段为0.1 Hz~n kHz不同^[3], 为了与常规音频大地电磁法区分, 定义为高频大地电磁法。目前, 高频大地电磁法可采用由GEOMETRICS和EMI公司联合开发生产的电磁观测系统EH-4电导率张量测量仪来实现。笔者近年来在我国东南部多个省市选择了近百个有代表性的测点, 实测了10 Hz~100 kHz频段电磁场, 计算出我国东南部天然电磁场平均振幅强度(图1), 结果显示在20 kHz~100 kHz频率范围内, 电场、磁场信号强度整体减弱, 电磁场信号相关度降低、重复性差。这是由于高频段信号易受人文噪声干扰, 如高压电线、电气铁路、电信接地设备等, 复杂多样的噪声影响往往使观测的重复性变差, 阻抗的估算分散, 不能客观地反映地下电性分布的情况, 甚至可能会得出错误的结果。

王书明等人^{[4, 5]}认为, 大地电磁信号是非高斯的、非最小相位和非线性的, 即是一个“ 三非” 信号, 经典的谱估计方法丢失了相位信息, 且对高频人文噪声的抑制不理想。笔者尝试采用高阶谱估计方法重构原始时间序列, 压制有色高斯噪声, 进而计算电阻率、相位, 从应用效果上看, 该方法在压制有色噪声和提取实际地质信息方面优于传统方法。

	Figure Option View Download New Window
	图1 中国东南部高频电磁场强度平均振幅

1 HMT信号的高阶谱的分析

设信号x(k)为零均值的k阶平稳随机过程, 则该过程的k阶矩m_kx(τ ₁, …, τ _k-₁)定义为^[6]

$\begin{matrix} \begin{matrix} m_{kx} (τ_{1}, \dots, τ_{k - 1}) = mon {x (n), x (n + τ_{1}), \dots \\ \dots, x (n + τ_{k - 1})}; \end{matrix} \end{matrix}$

该过程的k阶累积量 $\begin{matrix} C_{k}^{x} \end{matrix}$ (τ ₁, …, τ _k-₁)定义为

$\begin{matrix} \begin{matrix} C_{k}^{x} (τ_{1}, \dots, τ_{k - 1}) = E {x (n), x (n + τ_{1}), \dots \\ \dots, x (n + τ_{k - 1})} - E {g (n), g (n + τ_{1}), \dots \\ \dots, g (n + τ_{k - 1})}; \end{matrix} \end{matrix}$

其中:k≥ 3, g(n)是一个与x(n)具有相同功率谱的高斯过程。

由此, x(k)信号的三阶累积量为

$\begin{matrix} c_{3}^{x} (τ_{1}, τ_{2}) = E [x (k) x (k + τ_{1}) x (k + τ_{2})], \end{matrix}$ (1)

其中E{· }表示数学期望。经傅里叶变换, 式(1)变为

(2)

式(2)为x(k)的三阶谱^{[7, 8]}, 也称作双谱, 式中ω 为频率。三阶谱具有对称性, 即:

$\begin{matrix} c_{3}^{x} (ω_{1}, ω_{2}) = c_{3}^{x} (ω_{2}, ω_{1}) 。 \end{matrix}$ （3）

由式(2)可知双谱的相位谱为

$\begin{matrix} φ_{3}^{x} (ω_{1}, ω_{2}) = ϕ_{x} (ω_{1}) + ϕ_{x} (ω_{2}) - ϕ_{x} (ω_{1} + ω_{2}), \end{matrix}$ （4）

振幅谱为

$\begin{matrix} | C_{3}^{x} (ω_{1}, ω_{2}) | = | X (ω_{1}) | | X (ω_{2}) | | X (ω_{1} + ω_{2}) |, \end{matrix}$ （5）

其中ϕ _x(ω )和|X(ω )|分别为x(k)的相位谱和振幅谱。

2 高阶谱估计重构功率谱

2.1 幅值估计

由式(5)两边取对数可得^[5]

$\begin{matrix} \ln | X (ω_{1} + ω_{2}) | = \ln | C_{3}^{x} (ω_{1}, ω_{2}) | - \ln | X (ω_{1}) | - \ln | X (ω_{2}) | 。 \end{matrix}$ (6)

将式(6)作变量替换, 即令ω ₁+ω ₂=i, ω ₂=j, ω ₁=i-j, 可以得到离散表达式

$\begin{matrix} \ln | X (i) | = \ln | C_{3}^{x} (i - j, j) | - \ln | X (i - j) | - \ln | X (j) |; \end{matrix}$ （7）

当取i=0, j=0时, 可得

$\begin{matrix} \ln | X (0) | = \frac{\ln | C_{3}^{x} (0, 0) |}{3}; \end{matrix}$ （8）

当取i=1, j=0时, 可得

$\begin{matrix} \ln | X (1) | = \frac{\ln | C_{3}^{x} (1, 0) | - \ln | X (0) |}{2} 。 \end{matrix}$ （9）

由此, 可计算振幅谱|X(0)|、|X(1)|, 然后利用式(7)可以递推计算|X(n)|(n=1, 2, …, N)。

2.2 相位估计

双谱相位估计方法有多种^[5], 本文采用的是BMU算法。Brillinger最早提出了一种由φ (ω ₁, ω ₂)递推求ϕ (ω )的方法, 其出发点关系式为

$\begin{matrix} ϕ (ω) = \frac{1}{ω} [2 \int_{0}^{\infty} ϕ (λ) dλ - \int_{0}^{\infty} φ (λ, ω - λ) dλ] 。 \end{matrix}$ （10）

为了方便运算和提高计算精度, Matsuoka与Ulrych推导出了式(10)的离散形式, 在区域ω ₁+ω ₂=ω (0≤ ω ₁≤ ω ; 0≤ ω ₂≤ ω ), 进一步在ω ₁=[0, ω ]求和, 令ω ₂=ω -ω ₁可得

为了方便计算, 令Δ ω =1, ω ₁=i, ω ₂=j, ω =n, 可得

（11）

其中,

因此,

将上式改为递推形式, 定义

则有

其中, n=N对应于ω =π 。公式的初始值为

至此, 我们已经从大地电磁信号的三阶谱中分别提取出了振幅和相位, 得到信号的频谱, 经过反傅里叶变换重构信号, 进而进行功率谱估计求得各阻抗张量分量、视电阻率及相位信息。

3 理论与实测数据处理分析

3.1 数值模拟分析

采用数值模拟合成数据检验高阶谱技术压制HMT信号随机噪声的效果。数字模拟高频大地电磁高频信号:采样点数4 096, 采样频率12 kHz, 其中包含频率分别为30、50、100、200、300、1 000, 2 000 Hz的正弦信号, 利用matlab程序rpiid函数加有高斯噪声、非线性漂移和直流分量。

图2a为数值模拟的合成信号, 图2b~d分别为原始信号、重构信号傅里叶变换的振幅谱和相位谱。由图2b可知, 尽管原始信号经傅里叶变换能分辨到7个基频, 但是其抗干扰性很差, 100 Hz以上的频谱受噪声干扰严重, 其中1 000、2 000 Hz基本淹没在噪声中。由图2d可知, 原始信号的相位谱噪声的干扰严重, 基本无法有效地分辨各基频。由于高斯噪声的高阶统计量为0, 所以通过高阶谱重构的信号是压制高斯噪声的有效手段, 图2c、e是通过高阶谱重构信号的相应傅里叶变换振幅谱、相位谱, 能准确、清晰地分辨出设计的7个基频。

通过数值信号模拟可以证明, 高阶谱重构信号能够压制随机高斯噪声, 提高数据信噪比。

	Figure Option View Download New Window
	图2 数值模拟大地电磁信号

3.2 实测资料处理结果

大地电磁测深实际采集数据经常出现频点离散、误差棒较大等现象, 当前数据处理软件不能有效压制随机干扰噪声以及大地电磁信号非稳态特征是导致这一现象的原因之一, 通过数值模拟可知高阶谱重构信号技术在压制高斯噪声、提高信噪比方面有明显效果^[9]。图3为应用高阶谱信号重构技术进行大地电磁数据处理的基本思路。

	Figure Option View Download New Window
	图3 基于高阶谱重构技术数据处理流程

图4为实测数据通过高阶谱重构时间序列所得结果与传统处理结果对比。图4a为某隧道实测大地电磁信号, 电阻率曲线在10~100 Hz和2 MHz~10 MHz段出现锯齿状波动, 特别是在2 MHz~10 MHz段, 高、低阻交替出现, 整体看该点数据离散度大且误差棒大, 数据质量不可靠。图4b为高阶谱重构的大地电磁信号, 曲线整体趋势基本与图3a一致, 但其频点数据均匀平滑分布, 数据没有明显的跳跃, 误差也在允许范围之内, 特别是在2 M~10 MHZ频带范围内明显剔除了由随机高斯噪声引起的高阻, 数据质量明显优于原始信号。图4c为传统数据处理结果, 电阻率等值线较为凌乱, 浅地表和深部局部有“ 牛眼状” 高阻闭合圈出现, 基本无法有效提取地质特征。图4d为通过高阶谱重构的大地电磁信号, 等值线分布特征与图4c有一致性, 但其等值线平稳光滑分布, 很大程度上去除了浅地表及深部局部“ 牛眼状” 高阻闭合圈现象。图4d中, 电阻率断面在横向上分布较平稳, 与沉积岩宏观地质背景吻合, 纵向上整体呈高阻— 低阻— 高阻的趋势, 能有效地分辨出二叠系砂岩中的煤系地层(见图4d中虚线所示), 推测煤层顶板发育深度在地下170 m左右, 与后期钻孔揭示深度相吻合。

	Figure Option View Download New Window
	图4 实测原始信号和重构信号数据处理结果对比

4 结论与建议

高频大地电磁信号在20 kHz~100 kHz频率范围内, 电场、磁场信号强度整体减弱, 电磁场信号相关度降低、易受干扰、重复性差, 本文提出的基于高阶谱技术重构大地电磁信号, 通过数值合成数据和实测数据处理分析, 可以有效地压制视电阻率曲线中常出现的频点离散、误差棒较大等现象, 有效压制随机高斯噪声, 提高信噪比。

利用高阶谱技术重构大地电磁信号仍然存在一些问题需要进一步研究, 如高阶谱的实现方法、相位估计的初值选取都会不同程度影响最终的数据处理结果。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	王烨, 曹哲明, 汤井田, 等. 铁路隧道工程勘察中高频大地电磁测深应用效果研究[J]. 工程地质学报, 2005, 13(3): 424-428. [本文引用:1]
[2]	杨毅明, 化希瑞. 高频大地电磁测深在深埋隧道中应用实例分析[J]. 工程地球物理学报, 2009(S1): 37-41. [本文引用:1]
[3]	汤井田, 何继善. 可控源音频大地电磁测深法及其应用[M]. 长沙: 中南大学出版社, 2005. [本文引用:1]
[4]	王书明, 王家映. 大地电磁信号统计特征分析[J]. 地震学报, 2004, 26(6): 669-674. [本文引用:1]
[5]	王书明, 王家映. 利用高阶谱重构功率谱抑制高斯有色噪声[J]. 科学技术与工程, 2004, 4(2): 69-72. [本文引用:3]
[6]	张贤达. 时间序列分析-高阶统计量方法[M]. 北京: 清华大学出版社, 1996. [本文引用:1]
[7]	蔡剑华, 胡惟文, 任政勇, 等. 基于高阶统计量的大地电磁数据处理与仿真[J]. 中南大学学报: 自然科学版, 2010, 41(4): 1556-1560. [本文引用:1]
[8]	黄健英, 李录明, 罗省贤. 基于高阶谱的地震子波估计[J]. 成都理工大学学报, 2006, 33(2): 188-192. [本文引用:1]
[9]	化希瑞, 汤井田, 朱正国, 等. EH-4系统的数据二次处理技术及应用[J]. 地球物理学进展, 2008, 23(4): 1261-1268. [本文引用:1]

2005

0.0

王烨, 曹哲明, 汤井田, 等. 铁路隧道工程勘察中高频大地电磁测深应用效果研究[J]. 工程地质学报, 2005, 13(3): 424-428.

... 0 引言高频大地电磁法(HMT)是根据深埋隧道的工程地质勘探需要提出来的,工程勘察最关心的是近地表1 000 m内的地层地质情况^[1,2],要求大地电磁法观测频段为10 Hz~100 kHz,与传统的大地电磁法(MT)观测频段为0 ...

0.0

杨毅明, 化希瑞. 高频大地电磁测深在深埋隧道中应用实例分析[J]. 工程地球物理学报, 2009(S1): 37-41.

2005

0.0

... 1 Hz~n kHz不同^[3],为了与常规音频大地电磁法区分,定义为高频大地电磁法 ...

2004

0.0

王书明, 王家映. 大地电磁信号统计特征分析[J]. 地震学报, 2004, 26(6): 669-674.

... 王书明等人^[4,5]认为,大地电磁信号是非高斯的、非最小相位和非线性的,即是一个#cod#x0201c ...

2004

0.0

王书明, 王家映. 利用高阶谱重构功率谱抑制高斯有色噪声[J]. 科学技术与工程, 2004, 4(2): 69-72.

... 王书明等人^[4,5]认为,大地电磁信号是非高斯的、非最小相位和非线性的,即是一个#cod#x0201c ...

... 1 幅值估计由式(5)两边取对数可得^[5] ...

... 2 相位估计双谱相位估计方法有多种^[5],本文采用的是BMU算法 ...

1996

0.0

... _k-₁)定义为^[6] ...

2010

0.0

蔡剑华, 胡惟文, 任政勇, 等. 基于高阶统计量的大地电磁数据处理与仿真[J]. 中南大学学报: 自然科学版, 2010, 41(4): 1556-1560.

针对大地电磁(MT)数据具有非线性、非平稳的特征而不能满足传统功率谱方法要求的问题,从高阶统计量的理论入手,提出由信号的高阶谱恢复功率谱,再由恢复的功率谱估算MT响应函数的高阶统计量方法,建立通过高阶谱恢复功率谱的数学模型,并对合成数据和实测MT数据进行仿真分析.研究结果表明:利用高阶累积量对高斯噪声不敏感的特点,可实现在高斯白噪声和有色噪声污染下的瞬态信号频率与功率谱的正确估计;高阶统计方法消除了传统功率谱方法中大地电磁响应函数出现个别频点分散和形态异常的现象,使地质参数的估计精度和资料的可解释性得到明显提高.

... 式(2)为x(k)的三阶谱^[7,8],也称作双谱,式中#cod#x003c9 ...

2006

0.0

黄健英, 李录明, 罗省贤. 基于高阶谱的地震子波估计[J]. 成都理工大学学报, 2006, 33(2): 188-192.

摘　要：常规的地震子波估计是在假设子波为最小相位的基础上对地震记录进行处理。针对这一局限性，研究了基于地震记录高阶谱的地震子波估计方法。该方法无需对子波作任何假设，直接从地震记录的高阶谱中恢复地震子波的振幅谱和相位谱，进而重构地震子波。并对不同相位地震子波的理论模型和野外实际资料进行了处理，结果表明这一方法对于提高地震资料的分辨率具有可行性和应用潜力。

... 式(2)为x(k)的三阶谱^[7,8],也称作双谱,式中#cod#x003c9 ...

2008

0.0

化希瑞, 汤井田, 朱正国, 等. EH-4系统的数据二次处理技术及应用[J]. 地球物理学进展, 2008, 23(4): 1261-1268.

EH-4 system′s data processing mostly depends on IMAGEM which is manufacturer′s software. IMAGEM is inconvenient because of working in DOS interface, and has some disadvantages, such as low-resolution factor, asymmetrical frequency distribution, soit cannot effectively suppress disturbance. On the basis of studying EH-4 data processing, we ase a menthod to read time-sires file correctly, adopt periodogram means to recalculate auto-power spectrum, cross-power spectrum and tensor impedance component, and recalculate resistivity depend on tensor impedance component. We compile proprietary softwareEH-4 Data Processing Program. Because the quantity of reorganization impedance file is more 6 times than the original file, it enchances greatly vertical data density, resolution factor and credibility of data interpretation, According to actual EH-4 data analysis, the result is consistent with actual geological situation,.It proves that the program is feasible and valuable.

EH-4电导率成像系统的数据处理仍依赖于DOS界面下的数据采集和处理软件IMAGEM，不仅使用不便，而且有效频点少、频率分布不均匀、抗干扰能力差.本文在剖析EH-4系统数据处理流程的基础上，正确地读取了EH-4系统原始时间序列文件，采用周期图方法重新计算了自功率谱、互功率谱和张量阻抗分量，以及由阻抗定义的视电阻率，编制了专用软件EH4 Data Processing Program.由于重构阻抗文件的频点数约是原始阻抗文件的6倍，极大地增加了纵向数据密度，提高了分辨率和解释可靠性.实例应用表明了方法的正确、有效和可靠，具有推广价值.

... 2 实测资料处理结果大地电磁测深实际采集数据经常出现频点离散、误差棒较大等现象,当前数据处理软件不能有效压制随机干扰噪声以及大地电磁信号非稳态特征是导致这一现象的原因之一,通过数值模拟可知高阶谱重构信号技术在压制高斯噪声、提高信噪比方面有明显效果^[9] ...