跳到主要内容
Sandbox Physics
地震成像 · 边学边做实验

L4 · 第 8 章

读 Rawlinson 与 Spakman（2016）论地震层析的分辨率检验

请把 N. Rawlinson 和 W. Spakman 的论文 “On the use of sensitivity tests in seismic tomography” 打开放在旁边，发表于 Geophysical Journal International 第 205 卷，1221–1243 页（2016）。乌得勒支大学的机构库里有开放获取的全文。第 2 节用几行代数把合成检验和分辨率矩阵联系起来；第 3 节在同一个合成数据集上做了十四幅图的实验。读完后，你应该能对任何一篇论文里的分辨率检验说清楚：它取样的是分辨率矩阵的哪一部分，可能犯了哪一种常见错误。

Sandbox Physics ·

第 2.1、2.2 节：随时要对照的三个公式

式 6 是线性化方程组：数据等于矩阵乘模型，再加一个误差项；MIT 12.201 讲义第 4.19 节对直射线写出了同一个方程组，矩阵元素是各块内的路径长度。式 6 后面紧接着的那句话要读：误差项包含读数噪声、正演理论误差、线性化误差和参数化误差。后文的大部分论证，都在讨论某种检验能看到其中哪几项。

式 7 是代价函数，即数据失配加上由权衡参数加权的罚项；式 9 给出让它取最小值的广义逆。常见的误读是把分辨率当成射线本身的性质。式 9 表明广义逆取决于数据协方差和正则化，所以阻尼一变，式 13 的分辨率矩阵也跟着变。

式 13 把解拆成两项：分辨率矩阵乘真实模型，加上广义逆作用在误差上的部分。式 14 只用数据协方差就给出模型协方差。两项要分开看：分辨率说明真实结构怎样被模糊，协方差说明噪声怎样扩散。

d=Gm+ϵ  (6),G−g=(GTCd−1G+α2DTD)−1GTCd−1  (9),m^=R m+G−gϵ  (13)\mathbf d=\mathbf G\mathbf m+\boldsymbol\epsilon\ \ (6),\qquad \mathbf G^{-g}=\big(\mathbf G^{\mathsf T}\mathbf C_d^{-1}\mathbf G+\alpha^2\mathbf D^{\mathsf T}\mathbf D\big)^{-1}\mathbf G^{\mathsf T}\mathbf C_d^{-1}\ \ (9),\qquad \hat{\mathbf m}=\mathbf R\,\mathbf m+\mathbf G^{-g}\boldsymbol\epsilon\ \ (13)

第 2.3 节：为什么合成检验取样的是分辨率矩阵

式 15 和 16 是全文的核心。用同一个矩阵从检验模型生成合成数据，再用同一个广义逆反演，输出就是同一个分辨率矩阵乘检验模型，加上广义逆作用在合成噪声上的部分。

无噪声检验因此只测分辨率；只反演噪声的检验测误差传播；带噪声的检验把两者混在一起。论文里出现带噪声的棋盘图时，先弄清楚它想回答的是哪个问题。

ds=Gms+ϵs  (15),m^s=R ms+G−gϵs  (16)\mathbf d_s=\mathbf G\mathbf m_s+\boldsymbol\epsilon_s\ \ (15),\qquad \hat{\mathbf m}_s=\mathbf R\,\mathbf m_s+\mathbf G^{-g}\boldsymbol\epsilon_s\ \ (16)

第 2.5 节：八条推论里要标出的四条

第 (ii) 条：落在零空间里的检验模型完全恢复不出来，落在行空间里的则能完全恢复。作者在这里引用了 Lévêque 等人（1993）的理想化实验：射线相同，细棋盘被完美恢复，粗棋盘却几乎全在零空间里。棋盘恢复得好，并不证明其他尺度也分辨得好；检验只能暴露分辨率的不足。

第 (iii) 条：单个尖峰给出分辨率矩阵的一行和一列，相隔足够远的多个尖峰一次给出许多列。第 (iv) 条：检验必须使用真实反演的那个矩阵；在检验模型里重新追踪射线，会得到另一个矩阵，检验的就是别的东西了。第 (viii) 条：把恢复出的模型本身当作检验输入，如果分辨率矩阵接近幂等，得到的几乎还是原模型，这种检验会掩盖假象。

阅读时注意一处：第 3.6 节把最后这条论证称为第 vii 条，但在列表里它是第 (viii) 条。第 (vii) 条讲的是式 19 的叠加原理。

图 1 至图 5：棋盘与尖峰

图 1 是后面所有实验的基础。(a) 是真实模型，高斯异常最大正负 20%；(b) 是穿过它的射线；(c) 是 20 次非线性迭代后恢复的模型，走时残差均方根从 18.4 s 降到 0.35 s；(d) 是穿过恢复模型的射线。西北部覆盖好，往东南逐渐变差，最后几乎只剩一个方向。图 2 到图 5 的检验都沿用 (d) 的射线。

图 2 是细棋盘，图 3 是大小和幅度相同、但分布稀疏的尖峰。把两幅图并排对照：棋盘让人以为东南部完全分辨不了；尖峰则显示那里的大多数异常部分恢复了，只是被拉长，长宽比可达 5 到 6 比 1，而且能看出拉长的方向。要避免的误读是：把棋盘图上一片空白的区域当成毫无信息。

图 5 把尖峰检验的结果和真正的分辨率矩阵放在一起比较。模型有 3844 个未知量，作者用 Cholesky 分解算出整个矩阵约需 12 分钟；尖峰检验约 8 秒，结果基本一致。

图 6 至图 14：非线性、尺度、噪声与配色

图 6 和图 7 不再沿用观测射线，而是在检验模型中重新追踪射线。初至射线会弯向高速区，覆盖随之改变，棋盘沿对角线的拖尾明显加重。这正是第 (iv) 条的直观体现。图 8 估计了非线性误差项，并显示异常幅度减半后它明显变小。

图 9 把异常放大到约四倍再做检验，恢复变好了：在这个设置下，大结构比小结构更容易分辨。图 10 把恢复模型重新输入，几乎原样返回，对应第 (viii) 条的警告。图 11 改变阻尼和光滑，说明合成检验必须使用真实反演的正则化。图 12 加入 1 s 和 5 s 的高斯噪声；不加正则化时，5 s 的情形会过度拟合，凭空造出结构。图 13 显示参数化无法表示输入模型时带来的误差。图 14 显示，仅仅换一种配色，对分辨能力的印象就会不同。

摘要最后给出六条建议：用稀疏尖峰代替密集棋盘；沿用由观测约束的模型对应的数据覆盖；检验多种尺度；对模仿恢复结果的检验模型格外小心；用接近实际的噪声水平检验；条件允许时评估参数化误差。

完整例题

加阻尼的四单元网格不是幂等的

第 (viii) 条取决于分辨率矩阵有多接近幂等，即作用两次是否等于作用一次。取第六章的四单元网格：两条行射线、两条列射线，每个单元内长 1 km，数据误差取单位值，阻尼权重取 1。

数据项矩阵的特征值是 4、2、2、0。加阻尼后，每个特征值变成它自己除以它自己加一，所以分辨率矩阵的特征值是 0.8、0.667、0.667、0，每个对角元都是 0.533。只有特征值全为 0 或 1 的矩阵才幂等，这个矩阵不是。

在左上单元放一个单位尖峰。恢复图像在该单元为 0.533，两个相邻单元各 0.200，对角单元为负 0.133。按图 10 的做法，把这幅图像当作新的检验输入，输出是 0.382、0.160、0.160 和负 0.062，并没有返回输入。论文图 10 中几乎完全恢复，反映的是作者的设置和正则化，他们称之为“至少准幂等”。上面加阻尼的四单元网格显示，换一种设置可以离这个性质很远。

R=(GTG+I)−1GTG:  eig=45,23,23,0;R e1=(0.533, 0.200, 0.200, −0.133),R2e1=(0.382, 0.160, 0.160, −0.062)\mathbf R=\big(\mathbf G^{\mathsf T}\mathbf G+\mathbf I\big)^{-1}\mathbf G^{\mathsf T}\mathbf G:\ \ \mathrm{eig}=\tfrac45,\tfrac23,\tfrac23,0;\qquad \mathbf R\,\mathbf e_1=(0.533,\,0.200,\,0.200,\,-0.133),\quad \mathbf R^2\mathbf e_1=(0.382,\,0.160,\,0.160,\,-0.062)

读下一篇层析论文时的检查清单

找出检验模型：棋盘、尖峰，还是恢复模型本身？用了几种尺度？看合成数据是否用了真实反演的射线和矩阵，以及相同的阻尼、光滑和参数化。记下检验是无噪声、只有噪声还是带噪声，噪声多大。看看配色以及它在哪里饱和。最后问一句：论文是否把分辨率和模型不确定度分开展示了？两者各需要自己的检验。

带着问题动手

先预测，再改变，最后解释

  1. 在层析实验中选“台站环绕四周”，把真实模型先设为“棋盘”，再设为“尖峰”。预测在射线覆盖的边缘，哪一种图样能更清楚地显示拖尾方向。
  2. 在“棋盘”和“尖峰”之间切换真实模型，每次都把绿色探查格拖到边缘的一个单元上。
  3. 用第 2.5 节第 (iii) 条以及论文的图 2、图 3 解释两者的差别。
反演射线延迟，探查分辨率

这些实验使用平坦、各向同性的弹性或声学岩石，波速是虚构但接近实际的数值，都不对应真实的观测。每个实验在模型说明里写明模型和适用范围，并可导出读数背后的数据。

试着自己解释

1. 某篇论文展示了一个加了 3% 高斯噪声的棋盘检验，并把恢复出来的区域称为“分辨良好”。这幅图混合了哪两种效应？怎样把它们分开？

展开答案与解释

按式 16，输出等于分辨率矩阵乘棋盘，加上广义逆作用在噪声上的部分。无噪声检验单独给出前者；只反演噪声则单独给出后者。

2. 为什么一个完美恢复的棋盘仍可能掩盖没有分辨的大尺度结构？

展开答案与解释

能否恢复，取决于模型的哪一部分落在矩阵的行空间里。Rawlinson 和 Spakman 引用的 Lévêque 等人（1993）的例子就是这样：射线相同，细棋盘完全恢复，粗棋盘却几乎全在零空间里。检验因此要覆盖多个尺度。

3. 算出完整的分辨率矩阵约需 12 分钟，尖峰检验约 8 秒。两者相差多少倍？为什么在图 5 中尖峰检验给出的信息几乎相同？

展开答案与解释

约 90 倍。相隔足够远的尖峰各自取出分辨率矩阵的一列，彼此几乎不重叠，所以一次反演就能得到许多列。

12⋅60 s8 s=90\frac{12\cdot 60\ \mathrm s}{8\ \mathrm s}=90

可以放在旁边读的来源

  1. MIT OpenCourseWare, 12.201 Essentials of Geophysics (Fall 2004), lecture notes, Chapter 4: Seismology

    第 4.6 节把弹性波动方程分成 P 波和 S 波，并给出两者的波速（式 4.34 和 4.37）。第 4.14 至 4.16 节由 Fermat 原理推出 Snell 定律（式 4.91 至 4.96），并说明射线参数就是走时曲线的斜率（式 4.98）。第 4.19 节提到，走时层析中矩阵的元素就是射线在各小块里的长度。

  2. N. Rawlinson and W. Spakman, On the use of sensitivity tests in seismic tomography, Geophysical Journal International 205, 1221–1243 (2016)

    第 2 节写出阻尼最小二乘解和分辨率矩阵（式 7 至 14），并证明无噪声的合成检验采样的是同一个分辨率矩阵（式 16）。第 2.5 节列出八条推论。第 3 节在一个有 3844 个未知数的面波例子上比较棋盘检验和尖峰检验，最后给出六条建议。

反馈这个页面有错误、看不明白，或者缺了什么？告诉我们
反馈类型

提交时会附带本页地址、Lab ID、网站语言和设备类型。请不要填写个人信息。反馈最多保存 12 个月。隐私政策