本文旨在记录处理台站数据中遇到的问题及其解决办法。

噪声成像原理

Weaver,Wapenaar等大牛从理论上推导出结论:在均匀散射波场中,任意两点间的
互相关与两点间格林函数仅存在振幅差异。

  • BL Weaver(2011) 阐述了相对简单且易于理解的证明过程 => Seismic Noise Correlation
  • Weaver & Lobkis(2001) 两种不同方法证明互相关函数是经验格林函数的良好估计 => Paper
  • Wapenaar & Fokkema(2006) 利用液态球体模型证明了上述结论 => Paper

噪声成像主要包含如下过程:

  1. 单台数据处理
  2. 噪声互相关与相关函数叠加
  3. 提取频散曲线
  4. 速度结构反演及误差评定

Tip:本文专注于记载及解决单台数据处理和函数互相关与相关函数叠加的内容

单台数据处理流程

单台数据指各个台站在研究时间段内的连续波形记录。

其处理流程包含如下步骤:

  • 基础处理及时间分割
  • 时间域归一化
  • 谱白化

基础处理及时间分割

该处理旨在从地震仪原始观测记录中还原出我们感兴趣的频段内的连续波形数据,并以此提高信噪比。

该过程包含如下步骤:

时间域归一化

原始观测记录中往往夹杂着诸多非噪声信号。地震仪记录中常包含源自地震,海洋的信号以及台站仪器自身的噪音信号等。这些信号将破坏台站周边噪声信号的随机性,使均匀散射场的前提无法得到满足,其中对数据影响最大的当属地震信号无疑。故,为尽可能减小地震对噪声相关结果的影响,我们需要在时间域对数据进行归一化。
Fig.1 地震活动性对噪声相关结果的影响
Fig.1描述了不同的地震活性的情况对各个台站对间互相关函数的影响。当台站完全满足均匀散射场要求时,台站间互相关函数应该关于0完美对称。但地震活动性的差异将使其偏向一侧。图片来自@Heiner Igel的课件。

时间域归一化目前发展出了如下方法:

  1. One-Bit方法
  2. 原始波形的滑动平均法
  3. 改进的滑动平均法
  4. Water-level方法
  5. 地震事件自动捕捉法
  6. 天花板方法

One-Bit方法

最先介绍One-Bit方法不是因为它是最好的归一化方法而是因为它是最简单粗暴的方法。该方法仅保留了相位信息,完全放弃了数据的振幅信息。其将所有正振幅令为1,负振幅令为-1,令振幅原本为零者不变。
该方法的实现原理写成伪代码即为:

1
2
3
4
5
6
if Amplitude > 0 :
Amplitude = 1
else if Amplitude == 0 :
Amplitude = 0
else if Amplitude < 0 :
Amplitude = -1

原始波形的滑动平均法

该方法利用某一个时间段内的振幅均值的倒数对该时间段中心点的振幅值进行加权处理。
其中涉及到三个步骤:

  • 遍历所有的数据点\(d_n\),确定某一数据点前后某个固定宽度(\(N\)个数据点)为归一化时间窗。即遍历选定时间\( n\)的振幅\( d_n \),确定时间以其为中心的\(2N+1\)个数据点为时间窗。
  • 对据时窗内的原始波形数据的绝对值取平均值。即$$ w_n=\frac{1}{2N+1} \sum\limits_{j=n-N}^{n+N} \vert{d_j}| $$
  • 该\(d_n\)比上\(w_n\)得到归一化后的数据。即:\(\tilde{d_n} = \frac{d_n}{w_n}\)

Tips:

  1. Bensen等发现该归一化时间窗长度为最大滤波周期的一半时滤波效果较好。
  2. 数据中的异常大或小的尖锐振幅将影响很大范围内的数据的归一化结果。

改进的滑动平均方法

该方法同样进行滑动平均,需要给定归一化时窗长度。但是,该方法利用滤波至地震频段后的波形而不是原始波形记录计算权重因子。
其中涉及四个步骤:

  • 将原始记录\(d_j\)滤波至地震频段记录\(\hat{d_j}\)。
  • 遍历所有的数据点\(d_n\),确定某一数据点前后某个固定宽度(\(N\)个数据点)为归一化时间窗。即遍历选定时间\( n\)的振幅\( \hat{d_n}\),确定以其为中心的\(2N+1\)个数据点为时间窗。
  • 对时窗内的滤波后数据的绝对值取平均值。即$$ \hat{w_n}=\frac{1}{2N+1} \sum\limits_{j=n-N}^{n+N} \vert{\hat{d_j}}| $$
  • 该\(d_n\)比上\(\hat{w_n}\)得到归一化后的数据。即:\(\tilde{d_n} = \frac{d_n}{\hat{w_n}}\)

Water-level方法

Water-level方法是一个迭代趋近算法。该算法规定了振幅阀值,当\(d_n\)超过阀值时才对\(d_n\)进行滑动平均。如此对数据多次迭代进行归一化,直至所有数据均小于该阀值为止。

地震事件自动捕捉法

自动捕捉方法与Water-level方法类似,规定了阀值。该方法认定大于阀值的振幅即为地震干扰波动,直接剪切掉该阀值后续30min的地震记录,使其均为零。

天花板方法

该方法同样规定了一个阀值,将所有大于该阀值的振幅均归一化为该阀值。

Notes:

  1. 上述最后三种方法中均涉及到阀值的选择,其有很大的不确定性。对于不同地区、不同台站的数据来说,该阀值均不一定一致。由于地震噪声成像往往涉及到多个台站的数据,故这些方法无法进行数据的自动处理。
  2. 原始波形的滑动平均同样存在一定一些缺陷。(见Tips

谱白化处理

真实的地震噪声中往往存在诸多源。其中某些源可能特别强,占主导成分,譬如主微震(约20~40s的周期)和次微震(约5~10s的周期)。
Fig.2 地震噪声中成分的频谱分析;该图为台站NM.SLM的PSD;
从Fig.2(来自IRIS)我们可以大致观察到,一般地震噪声的PSD 谱图并非如我们所愿各个成分的谱强度大致相等。于是为了得到白噪声,我们需要对时间域归一化后的地震噪声进行谱白化处理
Fig.3 谱白化处理原理图;
从Fig.3(来自xsgeo)我们便可大致明白了谱白化基本流程:

  • 对噪声数据做实函数的傅立叶变换=>numpy.fft.rfft
  • 对噪声谱域的数据进行滑动平均
  • 将谱域数据反变换至时间域得到谱白化后的信号=>numpy.fft.irfft

噪声互相关和相关函数叠加

根据噪声成像原理部分可知,我们进行了上述单台单一天次的数据处理后,需要对处理过的两个台站数据做互相关才能得到格林函数。为提高信噪比,我们需要做叠加。
该过程可表示如下:

1
2
3
4
5
correlation_function=[0,0,...,0] # 向量长度取决于互相关时的偏移时间
for day in stacking_day:
xcorr = Cross-Correlation(stationA,stationB) #台站A和B在day这一天做互相关
correlation_function += xcorr #相关函数叠加

修改历史

  • 2017-01-10: 初稿