【电子读书笔记】Bensen et al. 2007 用噪声获得面波频散

1. Summary

a. 现状：噪声层析快速发展

b. 介绍本文工作：介绍目前噪声数据处理发展到的现状，并通过例子讲解

i. 噪声数据处理四个阶段：

1) 单台站数据准备

2) 相关与时域叠加

3) 量频散曲线（通过对群速度与相速度做频率-时间分析）

4) 质量控制（误差分析、选合适的测量数据）

ii. 重要细节：

1) 可以用于许多种观测设置；全自动

2) 数据质量控制非常重要（识别并提出坏的测量数据；为好的数据给出质量评价）

3) 质控的主要标准：稳定性、稳健性（在数据测量环境的扰动水平下）

4) 可靠性的重要指标：可重复性（本文认为季节可重复性同测量误差的重要性相等）

5) Proxy曲线（联系起观测信噪比与平均测量误差）在没有长时间序列的时候，可以提供有效的测量误差期望

2. Introduction

a. 理论基础：散射场互相关可以提供对台站间格林函数的估计

b. 相关实践：

i. 上述理论对面波成立：通过使用长噪声序列或尾波地震记录在台站对之间的互相关

ii. 上述理论对地壳体波成立：使用随机噪声数据

iii. 案例介绍（Shapiro，最早的研究）→相关发展综述（简单列举）→指出问题（只用到了20s周期以下的微震频带）

c. 引出本文工作：

i. 本领域目前的亟待改善1：宽频带数据的使用

ii. 本领域目前的亟待改善2：处理流程不断改善，但却没有被很好地整理出来

iii. 本文工作：总结数据处理的现状（自从第一篇文章以来，Shapiro and Campillo 2004）

d. 随机噪声数据处理基本流程总结：

i. 单台站数据准备

ii. 相关与时域叠加

iii. 测量频散曲线

iv. 质量控制（误差分析、选择合适的测量等等）

→ 完成这些步骤后，会进行群速度、相速度层析、Vs模型反演（不在本文讨论范围内）

→ 本文方法专门用于Rayleigh波，但Love波的相关研究也正在出现（Cho et al. 2006）

→ 重点关注：灵活性与应用于不同数据集上的普适性（不同频带、大小测网、海陆地震计）

a) 海地地震计OBS：长周期频带被污染，由于水流和重力波作用下的海底形变带来的倾斜

b) OBS垂向分量该问题的缓解方法：用水平分量数据+一个用于差分的台站（当时还未见其应用到宽频噪声数据处理中）

e. 再次总结本文目的、意义。并总结后续各章节内容

3. Single Station Data Preparation

a. 总体介绍

i. 目的：准备单个台站的波形数据（需要强调噪声）

1) 去除地震信号

2) 去除仪器的instrumental irregularities 

ii. 主要流程：

1) 移除仪器响应

2) 去均值

3) 去趋势

4) 带通滤波

5) 时域归一化（可包含对波形的非线性改变）

6) 谱白化（可包含对波形的非线性改变）

iii. 特点：

1) 由于是单台站操作，所以相较于后续互相关、叠加来说，计算资源、耗时要少

b. 时域归一化（最重要的步骤）

i. 总体描述

1) 作用：削弱地震事件、仪器irregularities，台站附近非稳态（non-stationary）噪声源，对互相关结果影响

2) 重要细节：地震是自动化处理流程的最大障碍

a) 地震的发生不规则

b) 许多小地震没有被记录到地震目录中

c) 短周期面波的到时还不是很清楚

→ 结论：地震信号的去除必须是数据自适应的（data-adaptive），而不是按照目录（catalogue）来预先规定的（prescribe）

ii. 五种自动识别并移除地震信号等污染信号的方法：

1) One-bit normalization：正的设为1，负的设为-1 →提升信噪比；（等价于N=0，时窗长度为1的方法4）

2) Clipped waveform：限幅法，振幅的阈值设为给定天内信号的rms（均方根）

3) 截除法：如果振幅高于阈值，那么30分钟的波形被设为0（问题，阈值设置很随意，不同台站间还不一样）

4) 绝对值滑动平均归一化方法：计算固定长度时窗内的，波形绝对值的滑动平均，用该值的倒数为窗口中心的波形加权

a) N=0，相当于One-bit归一化

b) N→∞，相当于原始波形的振幅整体缩放了一下

c) 用带通滤波器最大周期的一半，做时窗长度最合适。而且这个长度可以在相当大的程度上变化，依然给出相似的结果。

d) 缺陷：不能外科手术般地精准移除狭窄的异常，因为该方法下，这种异常会不可避免地在较长时间内对权重起到影响。（One-bit无此问题）

5) 迭代Water-level normalization：为振幅设置一个“Water-Level”(是数据的日rms的若干倍)，超过Water Level的数据被减权（down-weighted）。迭代若干次，直至所有的波形都低于Water Level(6倍rms水平)。

a) 缺点：耗时

iii. 举例：原始数据+五种方法处理后的数据做互相关，比较结果

1) 原始数据、限幅法、截除法互相关结果：信噪比很差，很Noisy的互相关结果

2) One-bit、绝对值滑动平均归一化方法、迭代Water-level normalization：高信噪比的互相关结果

3) 介绍测试数据集（15 GSN stations）下，利用observed spectral SNR（第三节）方法比较的结果：

a) 结论：绝对值滑动平均归一化方法最好

iv. 为什么共推荐绝对值滑动平均归一化方法？

1) 对数据更好的灵活性和适应性

a) 举例：地震活动性高的区域，一些小地震振幅和噪声水平差不多，直接用绝对值滑动平均归一化方法，去不掉这些小地震（可以从低通滤波结果看出）

b) 灵活性和适应性：如果在低通滤波后的序列上能够看到明显的小震，那么就在低频的带通滤波后的序列上计算权值，应用到原始数据上

c) 结果：有小地震时，波形被显著压制（down-weight），对应着带通滤波结果中地震信号被有效移除，于是污染得到了改善。

v. 被时域归一化漏过的地震信号对互相关结果有什么影响？

1) 主要体现为假的先驱信号（例如一些时间靠前的高振福信号）

c. 谱归一化（或谱白化）

i. Earth Hum、Primary and Secondary Microseism peak：噪声谱中的一些特征（Rhie & Romanowicz 2004）

ii. 一些单频的噪声源：会domain（主导）互相关的结果（at positive lag）→可通过对互相关结果做谱分析发现是单频噪声源影响；可以做噪声源定位

iii. 谱归一化：可以较好地改善上述问题

→ 窄带滤波：也能改善该问题。但可能会对频散测量造成一些问题？

iv. 谱归一化也可以通过改善不同台站的频谱强度不一致的问题，从而使得频散测量效果更好

→ 频谱中长周期信号的振幅更高，因为它们在长距离的传播中更加相关

4. Cross-correlation, Stacking and signal Emergence

a. 问题：台站间距离与是否能够获得可靠的测量结果的问题（太短/太长的距离都不行？）

→ 策略：在所有的台站对之间都先进性相关，共有n(n-1)/2组

b. Cross-correlation方式：daily，在频率域

c. Stacking：

i. 含义：b.的结果变到时间域，然后一个个叠加起来（也可以在频率域叠加，从而不用省去变换回时间域的步骤）

→ 推荐方法：逐级叠加：daily→weekly→monthly→yearly

→ 存储时间序列：有正坐标和负坐标；存储时间序列长度决定于台站间距与群速度

a) 正lag部分：叫做causal signal（因果信号）

b) 负lag部分：叫做acausal signal（非因果信号）

c) ↑代表从台站间的相对的不同方向传播的波。

d) 性质：如果随机噪声源随方位均匀分布，那么因果信号和非因果信号是完全一样的。

e) 处理方法：两支变一支（averaging，取平均）

ii. 更长时间的叠加，可以提升信噪比（SNR ratio）（图示说明见图10）

iii. 信噪比的频率依赖性：是一个可以一定程度上表征出面波信号随着参与叠加的序列增多而emerge的量（→第四节：用于QC）

1) 介绍了一种：'spectral' SNR measurement

→ 引号的来源：实际上还是在时间域进行测量的

→ 测量方法：细节还有些没看懂，在第10页

2) 需要更多的研究来探索，哪一种SNR的定义方式对频散测量更有帮助

3) SNR关于频率与叠加时间长度、（台站间距）的规律的详细研究：

a) 指数律：SNR=At^(1/n)

b) 频率对应的噪声信号弱，则n大，SNR增长慢；在例如Earth-Hum等噪声强的频带，SNR增长更快

c) 这种指数律对噪声信号成像非常有用（见质控节）。但关于这种指数律本身的behavior，以及它的geographic variability，有待于进一步研究。

5. Dispersion Measurement

a. 经过步骤4.叠加完成后：结果是EGF（经验格林函数）→可以测量群速度和相速度（方法：时间频率分析FTAN，frequency-time analysis）

b. 分析方法介绍（FTAN，公式见第11页）

i. 傅里叶变换：注意变换时乘的是正指数（跟数学教材上不同，为什么？）

ii. analytic signal+高斯窄带滤波（滤波参数α可调，其效果跟时频域分辨率有关；通常根据距离来设定。）

iii. 反傅里叶变换→振幅谱（encelope function）、相位谱（Phase function）

iv. 群速度到时，根据振幅谱峰值给出

v. Instantaneous frequency的定义：根据相位谱时间变化率给出→用来替换ω0，主要是出于工程实际的考虑

c. Group Velocity Measurement：包括八个步骤

i. 第一步，获得Frequency-time image（FTAN image）

1) 公式：log⁡〖|A(t,ω_0 )|^2 〗

2) 群速度代替时间，周期代替速度

ii. 第二步，Track Dispersion ridge →（获得） Raw Group Speed Curve

iii. 第三步，phase-matched filtering method，去除波形中可能的污染

iv. 第四步，applying anti-dispersion filter to get undispersed signal

v. 第五步，识别污染信号，并且将其从undispersed signal中移除

vi. 第六步，重建提出污染信号后的信号

vii. 第七步，算cleaned信号的FTAN图

viii. 第八步，在cleaned FTAN图track频散ridge

→ 与前人方法的比对：

1) 前人：分析人员手动做步骤3-5

2) 改进：自动化处理步骤3和步骤5

a) 步骤3：自动化处理的问题 → 噪声数据避免了一些问题（多路径效应，较少的谱中holes），方法改进解决了一个问题（jump）

b) 步骤5：介绍了这一步处理的目的——去掉了一些假的前兆信号，使得FTAN图更平滑

→ 介绍Phase-matched filtering（步骤3-8）的缺陷：不适用于窄带的波形处理。因此只推荐在宽频带信号的处理中使用

d. Phase speed measurements

i. 背景：

1) 分析envelope function，可以得到群速度曲线

2) 相速度不能从群速度直接得到

3) 群速度可以通过相速度计算得到

ii. 群速度与相速度：

1) 群速度∂ω/∂k

2) 相速度ω/k

iii. 引入相位观测信息、群速度到时的观测信息→eq.10，相速度曲线

1) 该方程使用受到相位不确定性以及初始相位的制约

a) 初始相位

i) 真实地震：初始相位来自于CMT解（CMT：Centroid Moment Tensor，质心矩张量？）

ii) 噪声方法：初始相位认为是0

b) 不确定性

i) 不论什么情况，都有2πN的不确定性

→ 可以通过使用全球3D模型，或相速度图来找到N （例如，长周期T>40s，可能定N比较准。也可能要求T>100s比较准）

→ 第5节给出，频散测量周期最多使用到Δ/12，于是T>40s要求台站间距500Km

→ 噪声方法：短周期、短距离，因此相位不确定性可能难以通过这些直接的方法求得

ii) 噪声互相关中的其他不确定性来源

→ π/4项：随方位的均匀弥散场会有的不确定性。如果不均匀，shift还会不一样+有频率依赖？

→ Lin的文章：用垂向，大概还是π/4

→ 总之还需要更多的研究

6. Quality Control

a. 背景：

i. 互相关数据的量随台站数平方增加→要尽可能减少人工操作→人工操作少，自动质控问题就更重要

ii. 计算数据质量的相关指标，并且确定容忍标准

b. 初步质控标准：

i. 可靠频散测量要求：∆>3λ=3cτ

1) ∆震中距、λ波长、τ周期

2) 由于相速度c~4km/s，因此τ<∆/12

3) 对密集台阵而言，是一个较为严格的约束

4) 中长周期信号的测量：永久站点（基岩上建立）的重要性

5) 考虑到目前对于基于噪声互相关的相速度测量的经验较少，因此该准则可能还会被放宽

ii. 在上一步约束之外，还需要确定频散测量的质量如何

1) 角度1：与真实数据比较。需要地震刚好发生近似在台站下方（Fig. 14, 15，coincidence！）

a) 基阶瑞利波到时基本一致

b) （经过一堆复杂的处理，总之可以论证从观测上来说，噪声互相关信号得到的东西，跟真实地震信号是差不多的，'confidence'）

2) 角度2：数据的稳定性

a) 数据在扰动下的稳健性

b) 空间集群/时间重复的观测对于测试a)很有用（→可以给出标准差？之类的可用于表征误差的统计量）

i) 空间集群（clustering）：

• 例如，发生位置差不多的一系列地震→可以用来估计地震频散测量的不确定性

• 类似地，在噪声互相关中，高密台阵有空间集群的效果（Fig. 16）

> 集群分析：台阵集群，小张角

> 目前的作用：长路径测量的平均不确定性估计；某些测量子集的异常数据剔除

ii) 时间重复（repeating）：

• 对于噪声互相关方法而言，季节性（即有repeatable的意思）噪声源变化可以帮助估计测量不确定度（Fig. 17）

• 绘制一年内不同时间的互相关结果（要求：信噪比SNR大于某个阈值，例如10），结果的差异可以给出对不确定度的估计。

c) 对于不确定度无法计算的数据，予以剔除

i) 影响：信噪比差的情形下，可能要用更长时间的观测数据叠加（例如2年而不是1年）才可靠（叠加时间越长，结果越稳定；更长的叠加序列的话，做temporal repeating分析的时候，取的小时间片段可能也可以更长）

3) 角度3：信噪比SNR （proxy curve），也可以认为是角度2的一个延伸

a) 有时候不同台站之间同时运行的时间（Temporal Overlap）不够，就难以按前面的方法给出不确定性的估计 → 用信噪比SNR度量（Fig. 18）

b) 依据：（数据上→）SNR和不确定度之间的相关性分析（Fig. 18）；（理论上→）SNR跟叠加时长有关（Fig.12），叠加时长又跟测量不确定度的期望有关

i) 某数据集的分析结果：（结论：SNR>10时，SNR可能可以作为频散测量不确定度的不错的指标）

> spectral SNR和temporal repeating分析得到的不确定度（平均标准差）在SNR>10时有比较好的线性关系

> 在SNR<10时呈现强的非线性关系

> Fig.12：信噪比的幂衰减率 → 并不需要无穷的叠加，只要SNR高到能够给出一个合理的频散曲线测量效果就行了 → 对应某个最小叠加长度阈值L0

> 如果有大于L0的测量时间序列，更多的是可以用作temporal repeating分析，这比使用SNR proxy curves更好。

4) 阶段总结：more illustrative than definitive，需要做更多工作来guide their use

5) 角度4：测量结果可以对得上其反演的光滑层析结果 → cohere as a whole; agree with other accepted measurements （Fig. 19 and Yang 2007）

7. Summary and Conclusions

a. 现状：从噪声数据中提取估计的瑞利面波经验格林函数的流程，如今已经较为成熟，未来几乎不会出现根本性的modification。

b. 点出流程特点，总结四个步骤。指出自动化流程的必要性。

i. 单台站处理

ii. 互相关与叠加

iii. 频散测量

iv. 质控

c. 简要总结各个步骤要点

i. 单台站处理：主要做时域归一化（改善地震、仪器不规则相应、风暴等非稳态噪声对数据造成的污染）；推荐滑动绝对平均方法；推荐做谱白化

ii. 互相关与叠加：长时间数据→信噪比好；SNR的指数律；指数律中的指数大小与频段关系→30s-60s频段较为challenging

iii. 频散测量：FTAN；增加phase-matched filtering→帮助改善存在相邻的污染信号的情形，不能用于窄带波形；准则SNR>10，站间距大于三倍波长；相位确定问题

iv. 质控：首要准则→稳健性，时间可重复性，将季节可重复性同测量不确定性等同；短测量序列用SNR做质控的介绍

d. 展望：噪声层析方法历史不足三年，还有很多重要工作需要做

i. 为了提取更好的相速度，与源分布相关的相位不确定性

ii. 勒夫波

iii. 代理曲线proxy curvve