try ai
科普
编辑
分享
反馈
  • 走时层析成像

走时层析成像

SciencePedia玻尔百科
核心要点
  • 走时层析成像通过测量地震波的走时来创建地球内部的图像,地震波在高温、低密度的物质中传播较慢,而在低温、高密度的物质中传播较快。
  • 该技术在数学上被表述为一个大型反演问题 (d=Gm\mathbf{d} = G\mathbf{m}d=Gm),旨在寻找能够最好地解释观测到的走时数据 (d\mathbf{d}d) 的地球模型 (m\mathbf{m}m)。
  • 由于射线覆盖不均匀以及波传播的平滑特性,层析成像是一个不适定问题,需要使用正则化技术来获得稳定且地质上合理的解。
  • 层析成像图像的质量和局限性通过模型分辨率矩阵进行评估,该矩阵揭示了地质特征在反演过程中是如何被模糊或涂抹的。
  • 层析成像的底层数学框架是通用的,与机器人学中的同步定位与建图 (SLAM) 等其他领域的方法有着深刻的联系。

引言

我们如何能看见一个无法打开的行星的内部?这个问题驱动了地球物理学家数十年,并且是走时层析成像的核心。就像通过敲击一个密封的盒子来猜测其内容物一样,地震学家利用地震产生的“敲击”和全球传感器记录到的“回声”来构建地球深部内部的图像。其核心思想很简单:地震波从震源传播到接收器所需的时间揭示了它所穿过的物质的性质。然而,将这些时间测量值转化为我们星球的连贯图像,是一个复杂的科学和数学难题。本文旨在弥合简单概念与其复杂执行之间的知识鸿沟,全面概述了走时层析成像的工作原理,从第一性原理到高级应用。第一部分“原理与机制”将解析其数学基础,解释我们如何将走时转化为一个庞大的方程组,并应对解决这个不适定反演问题的挑战。随后的“应用与跨学科联系”部分将探讨如何使用这些模型来揭示地球内部的动力学过程,科学家如何处理现实世界中不完美的数据,以及这些强大的思想如何与机器人学等看似遥远的领域联系起来。

原理与机制

想象一下,你有一个密封的木箱。你不允许打开它,但你想知道里面有什么。你会怎么做?你可能会在不同的地方敲击它并仔细聆听。空洞的回声可能暗示着一个空的空间,而沉闷的响声则可能表明有固体物质。通过从各个侧面收集足够多的这类“敲击-聆听”测量,你就可以开始绘制出箱子内部的粗略地图。

走时层析成像就是这个简单想法在行星尺度上的放大版。“敲击”是地震或受控爆炸,而“聆听者”是散布在全球各地的数千个地震传感器(地震仪)。我们所聆听的信息是走时:地震波从震源传播到接收器需要多长时间。其核心原理异常简单:波在更致密、更刚性的材料中传播得更快,而在密度较低、温度较高或部分熔融的材料中传播得更慢。通过精确测量这些走时,我们就可以开始绘制出地球深部的快速和慢速区域,从而有效地为我们的星球进行一次CAT扫描。

从平滑地球到块状世界

将走时与地球结构联系起来的基本物理定律非常优美。波沿路径 γ\gammaγ 传播所需的时间 TTT 是材料​​慢度​​沿该路径的积分:

T=∫γs(x) dℓT = \int_{\gamma} s(\mathbf{x}) \, d\ellT=∫γ​s(x)dℓ

这里,s(x)s(\mathbf{x})s(x) 是空间中每一点 x\mathbf{x}x 的慢度,它就是波速的倒数,s=1/vs = 1/vs=1/v。将慢度在整个路径长度 dℓd\elldℓ 上进行积分(或求和),就得到了总走时。这个积分方程是我们工作的“罗塞塔石碑”;它将我们能够测量的量 (TTT) 与我们想要知道的量 (s(x)s(\mathbf{x})s(x)) 联系起来。

当然,地球是一个连续、复杂的物体,我们不可能确定每一点的慢度。为了让计算机能够处理这个问题,我们必须简化或​​离散化​​我们的模型。想象一下,在地球内部覆盖一个巨大的三维网格,将其划分为数百万个小方块,或称“体素”。然后我们做一个简化假设:在每个方块内,慢度是恒定的。这被称为​​分段常数参数化​​。

通过这种简化,我们那个优美但困难的积分方程就变成了一个简单的求和。单个射线的总走时现在是它在所穿过的每个方块中花费的时间之和。在方块 jjj 中花费的时间就是射线在该方块内的路径长度 LijL_{ij}Lij​ 乘以该方块的慢度 sjs_jsj​。对所有方块求和,射线 iii 的走时为:

Ti=Li1s1+Li2s2+⋯+Linsn=∑j=1nLijsjT_i = L_{i1}s_1 + L_{i2}s_2 + \dots + L_{in}s_n = \sum_{j=1}^{n} L_{ij}s_jTi​=Li1​s1​+Li2​s2​+⋯+Lin​sn​=j=1∑n​Lij​sj​

如果我们有 mmm 个这样的射线路径测量值和模型中的 nnn 个方块,我们就可以将其写成一个巨大的线性方程组。这就诞生了位于层析成像核心的著名矩阵方程:

d=Gm\mathbf{d} = G\mathbf{m}d=Gm

我们不要被这些符号吓到。它代表着非常具体的东西:

  • d\mathbf{d}d 是我们的​​数据向量​​,一个包含我们测量到的所有走时的长列表。
  • m\mathbf{m}m 是我们的​​模型向量​​,另一个包含我们想要发现的每个方块的所有未知慢度值的长列表。
  • GGG 是​​设计矩阵​​或​​雅可比矩阵​​。它并不神秘;它只是一个巨大的数字表格,其中每个条目 GijG_{ij}Gij​ 是射线 iii 在方块 jjj 中的路径长度 LijL_{ij}Lij​。它是将我们的地球模型 (m\mathbf{m}m) 与我们期望看到的数据 (d\mathbf{d}d) 联系起来的几何映射。

侦探的困境:解开谜题

我们现在面临一个经典的侦探问题。我们有证据 (d\mathbf{d}d) 和一套将证据与嫌疑人联系起来的规则 (GGG)。我们需要找出罪魁祸首——即真实的慢度模型 (m\mathbf{m}m)。这就是​​反演问题​​的本质。

在理想世界中,我们可以直接对矩阵 GGG 求逆来找到 m=G−1d\mathbf{m} = G^{-1}\mathbf{d}m=G−1d。但现实世界是混乱的。我们的测量从不完美;总有来自仪器误差、背景振动以及模型简化的噪声。所以真实的关系更接近于 d=Gm+ε\mathbf{d} = G\mathbf{m} + \boldsymbol{\varepsilon}d=Gm+ε,其中 ε\boldsymbol{\varepsilon}ε 是一个随机误差向量。

既然我们无法找到一个完美的解,我们就必须寻求能够最佳解释我们数据的模型 m\mathbf{m}m。最自然的方法是​​最小二乘法​​。我们寻找一个模型 m\mathbf{m}m,它能最小化我们的预测 GmG\mathbf{m}Gm 与实际测量值 d\mathbf{d}d 之间的差异——即“残差”。具体来说,我们最小化误差的平方和,这个量由平方范数 ∥Gm−d∥22\lVert G\mathbf{m} - \mathbf{d} \rVert_{2}^{2}∥Gm−d∥22​ 给出。

这个选择并非任意。事实上,它有深刻的理由。如果我们假设误差 ε\boldsymbol{\varepsilon}ε 是独立的并服从高斯分布(经典的“钟形曲线”),这对于许多小的、不相关的误差源来说是一个合理的假设,那么最小二乘解也就是​​最大似然估计​​。它是最有可能产生我们观测到数据的那组模型。为了使该方法产生一个单一、唯一的答案,一个数学条件是我们的矩阵 GGG 必须具有满列秩,这意味着我们的模型参数中没有冗余。

反演的陷阱:为何层析成像是困难的

所以,我们有了一个计划:建立矩阵方程并找到最小二乘解。会有什么问题呢?事实证明,问题很多。当我们实际尝试这样做时,得到的地球模型常常看起来像一团充满噪声、毫无意义的混乱。这是因为层析成像是​​不适定问题​​的一个经典例子,也正是在这里,该领域的真正挑战与魅力得以显现。

平滑的诅咒

第一个罪魁祸首是我们正演模型本身的性质:T=∫s dℓT = \int s \, d\ellT=∫sdℓ。积分是一个平均过程,它会使事物平滑化。想象一个慢度场,它在高值和低值之间快速振荡,就像一个棋盘格。穿过这个场的射线将会平均掉快速和慢速的部分。最终的走时可能与穿过一个完全均匀介质的射线的走时无法区分。换句话说,我们的正演模型天生就会滤除掉精细的细节。而反演问题必须做相反的事情:它必须对数据进行“去平滑”或微分,以恢复这些细节。这本质上是一个不稳定的操作。就像试图把炒好的鸡蛋复原一样,它会极大地放大我们数据中存在的任何微小噪声。

我们视野中的空白

第二个,也许更直观的问题是,我们的视野是有限的。地震和地震仪并非无处不在。这导致了​​不完全和不均匀的射线覆盖​​。地球地幔的某些区域被成千上万条射线纵横交错地穿过,而其他区域,如南半球深处,则是地质盲点。

为了以最鲜明的形式看到其后果,让我们想象一个微小的、简化的 2×22 \times 22×2 世界,其中有四个我们想要绘制的慢度块。假设我们的实验是有限的:我们只能纯粹水平和纯粹垂直地发射射线。我们测量每行和每列的慢度之和。 现在,考虑在慢度模型中添加一个“棋盘格”模式:我们将一个很小的量 δ\deltaδ 加到左上角 (s11s_{11}s11​) 和右下角 (s22s_{22}s22​) 的单元格,并从另外两个单元格 (s12s_{12}s12​ 和 s21s_{21}s21​) 中减去相同的量。我们的测量值会发生什么变化?让我们检查顶行:其总走时的变化是 (δ)+(−δ)=0(\delta) + (-\delta) = 0(δ)+(−δ)=0。没有变化!左列的变化是 (δ)+(−δ)=0(\delta) + (-\delta) = 0(δ)+(−δ)=0。同样没有变化!这个棋盘格模式,作为模型的一个真实特征,对我们的实验来说是完全不可见的。

这个不可见的模式是位于我们矩阵 GGG 的​​零空间​​中的一个向量。它的存在意味着解不是唯一的,而是存在一个无穷的解族(我们原始的模型加上任意数量的棋盘格模式),它们都能完美地拟合数据。矩阵 GGG 是奇异的,其条件数为无穷大。这反映在法方程矩阵 GTGG^T GGTG 的特征值中;零空间的存在保证了至少有一个特征值恰好为零。

在真实的、大规模的地球中,情况很少如此完美奇异。但原理是相通的。总会存在一些慢度变化的模式——通常是小尺度的复杂结构——我们特定的射线网络对它们不是很敏感。这些对应于一个​​近零空间​​。用线性代数的语言来说,这些约束不佳的模式与矩阵 GGG 的非常小的​​奇异值​​相关联。解决反演问题实际上需要除以这些极小的数。数据中恰好与这些方向一致的任何噪声都会被放大,从而污染解。这正是问题​​病态​​的含义,或者说矩阵 GTGG^T GGTG 是​​近乎秩亏​​的。

物理学家的工具箱:从直线到弯曲的世界

不适定性的挑战是深刻的,但物理学家和数学家已经开发出一套复杂的工具箱来应对它们。这包括使我们的模型更真实,我们的方法更稳健。

世界不是直线的

我们开始时假设地震波沿直线传播。但真的是这样吗?作为光学基石的费马原理指出,波总是会沿着走时最短的路径传播。在均匀介质中,这是一条直线。但在速度变化的地球中,最快的路径通常是一条曲线,它会弯曲以远离慢速区域并朝向快速区域。例如,在一个慢度随深度增加的分层介质中,以一定角度发射的射线在向地表弯曲时会遵循一条优美的圆弧。

这意味着射线路径——即编码在我们矩阵 GGG 中的几何结构本身——依赖于我们试图寻找的慢度模型 m\mathbf{m}m!这种反馈循环使得问题从根本上是​​非线性​​的。简单的线性方程 d=Gm\mathbf{d} = G\mathbf{m}d=Gm 不再是故事的全部。如果慢度变化很小,射线弯曲可以忽略不计,我们可以采用线性的、直线射线的近似。但为了获得高保真图像,我们必须接受非线性。这通常通过迭代方法来完成:猜测一个模型,追踪穿过它的弯曲射线,根据数据失配更新模型,重新追踪射线,并重复此过程,直到过程收敛到一个自洽的解。

选择你的现实:慢度与速度

一个微妙但强大的选择是我们模型的​​参数化​​。我们一直使用慢度 (sss),因为它使得基本的正演问题是线性的:T=∑LjsjT = \sum L_j s_jT=∑Lj​sj​。如果我们选择对速度 (vvv) 建模会怎样?由于 s=1/vs=1/vs=1/v,走时将是 T=∑Lj/vjT = \sum L_j/v_jT=∑Lj​/vj​。突然之间,即使我们假设射线是固定的直线,问题从一开始就是非线性的!最初选择使用慢度是一个深思熟虑的决定,目的是使数学尽可能简单。这揭示了建模中的一个深刻真理:你选择如何描述世界,会极大地改变解决这个难题的难度。

初至为佳

在一个复杂、不均匀的地球中,波会以复杂的方式反射和折射,从震源到接收器可能有多条不同的路径,从而在地震图上产生多个不同的到时。我们在测量中应该使用哪一个?答案几乎总是​​初至​​。这并非为了方便;这是费马原理的一个深刻推论。虽然所有有效的射线路径都是“稳态”时间的路径,但只有第一个到时才能保证是​​全局最小​​走时路径。后续的到时,可能是从结构上反弹或穿过称为​​焦散​​的强聚焦区域而来,对应于走时泛函中的局部极小值或鞍点。通过只关注初至,我们将数据建立在一个稳健、单值且数学上稳定的原理之上。它将一个潜在的、混乱的、多值的问题转化为一个行为良好(尽管仍是不适定的!)的反演问题,其导数是稳定且定义明确的。

构建更好的网格

最后,即使是我们的离散化方法也可以改进。我们可以不在粗糙的、允许之间存在尖锐、不自然跳跃的恒定慢度块上定义慢度,而是在节点网格上定义慢度,并在其间的空间中进行平滑插值。这从一开始就为我们的模型强制施加了一定程度的光滑性,起到了一种内置物理直觉的作用。它含蓄地惩罚了那些通常是噪声产物的狂野、高频棋盘格模式,从而改善了问题的条件,并导向更稳定、地质上更合理的图像。

走时层析成像的历程是一个精彩的科学探案故事。它始于对回声计时的简单想法,带领我们穿越线性代数的优美数学、反演理论的深刻挑战以及波传播的深层物理原理。它证明了,只要我们谨慎而巧妙,仅凭地震图上几条微弱、晃动的线条,就能构建出一个对未知世界惊人清晰的图像。

应用与跨学科联系

至此,在我们的探索之旅中,我们已经揭示了走时层析成像的基本原理。我们已经看到如何构建问题,将一组时间测量值转化为一个方程组,并且我们已经探索了解决它所需的数学机制。但科学并非在真空中完成。这些思想的真正美妙之处不在于其抽象的优雅,而在于它们让我们能够做什么。现在,我们走出教室,走向世界,看看层析成像是如何在实践中工作的,应对其现实世界的局限性,并发现它与其他科学和工程领域的惊人联系。

洞见未见之艺:为地球成像

地震层析成像最主要也是最令人惊叹的应用是创建地球内部的地图。对于地球物理学家来说,它是窥探深部地幔的宏伟望远镜,那是一个在我们脚下数千公里、永远无法直接观测的领域。通过收集全球地震台站记录的数千次地震的走时,我们可以构建出地球内部结构的图像。

在这些层析成像图像中,我们看到的不是岩石和矿物,而是地震波速的变化。波传播速度比平均速度更快或更慢的区域,会显示为明显的异常。慢速异常可能指示一个比周围环境更热的区域,例如,一个正在上升的地幔热柱,它为夏威夷或冰岛等地的火山提供物质。另一方面,快速异常通常指向更冷、更致密的物质,例如在俯冲带俯冲回地幔的巨大大洋板块。这些图像不是静态照片;它们是一个动态的、对流的行星的快照,为驱动板块构造的引擎提供了我们所拥有的最直接的证据。一个精心设计的反演,即便是合成算例,也能展示该技术解析此类特征的能力,但同时也揭示了其挑战。

当然,现实世界是混乱的。我们的数据从不完美,一个稳健的科学方法必须坦诚面对其不完美之处。层析成像完美地展示了如何处理两种常见的数据缺陷:不确定性和离群值。

想象一下你在为一场比赛计时。对于一些选手,你可能有一个非常精确的秒表;而对于另一些选手,你可能需要瞥一眼远处的时钟,导致测量结果不那么确定。对这两种测量给予同等的置信度是愚蠢的。对于地震数据也是如此。地震波的到时有时可以被高精度地拾取,而有时信号充满噪声,拾取结果不确定。加权最小二乘法提供了解决方案。我们不是最小化简单的误差平方和,而是最小化一个加权和,其中每个误差的贡献都根据我们对该测量的置信度进行缩放。高置信度的测量获得大权重,将解拉向它,而低置信度的测量则被赋予较小的影响。这不仅仅是一个聪明的技巧;可以证明,在关于误差性质的合理假设下,该过程会产生*最大似然估计*——即在给定我们的数据及其不确定性知识的情况下,最有可能的模型。

但如果你的一位计时员犯了一个巨大的错误——一个真正的离群值,该怎么办?标准的最小二乘法会狂热地试图减小误差的平方和,从而会受到严重干扰。一个大的误差,平方后会变得巨大,算法可能会扭曲整个模型,只为了减小那一个坏数据点的影响。我们需要一种更稳健的方法,一种对极端错误不那么敏感的方法。这就是诸如迭代重加权最小二乘法 (IRLS) 等方法发挥作用的地方,通常使用像Huber惩罚这样的“惩罚函数”。其思想非常直观:对于小的、合理的误差,该方法的作用与标准最小二乘法完全相同。但当遇到一个非常大的残差——一个可能的离群值时——它实际上会说:“这个数据点可能是错的,我不会那么努力去拟合它。”它系统地降低这些离群值的权重,防止它们破坏整个图像。这是一种直接内建于数学中的自动化科学怀疑主义。

地图并非疆域:理解分辨率

一幅层析成像图像是一项了不起的成就,但理解其局限性至关重要。地图并非疆域。它永远是现实的一个简化、平滑且不完美的表示。我们必须始终问的核心问题是:“我们真正能看到什么?”

答案在于一个称为​​模型分辨率矩阵​​的概念,我们可以用 R\mathbf{R}R 来表示。在理想世界中,如果真实的地球由一组参数 mtrue\mathbf{m}_{\text{true}}mtrue​ 描述,我们的反演将得到 m^=mtrue\hat{\mathbf{m}} = \mathbf{m}_{\text{true}}m^=mtrue​。分辨率矩阵将是单位矩阵。实际上,我们估计的模型是真实情况的一个模糊版本:m^=Rmtrue\hat{\mathbf{m}} = \mathbf{R} \mathbf{m}_{\text{true}}m^=Rmtrue​。分辨率矩阵 R\mathbf{R}R 就像我们层析成像“相机”的镜头,而它很少是一个完美的镜头。它会模糊和扭曲真实的图像。

为了理解这种模糊的性质,我们可以做一个思想实验。如果真实的地球是完全均匀的,除了在某一点上有一个微小的异常,它的图像会是什么样子?答案由分辨率矩阵的相应列给出。这一列就是​​点扩散函数​​ (PSF)。它在层析成像中相当于天文学家拍摄的单个遥远恒星的图像。一个完美的望远镜会显示一个光点;而一个真实的望远镜则显示一个小的、模糊的圆盘,可能周围还有光环。层析成像中的PSF精确地告诉我们,一个位置的异常在我们的最终图像中是如何被“涂抹”开并泄漏到邻近区域的。

分析这些PSF是评估层析模型质量最诚实的方法。几十年来,一种常见的做法是“棋盘格测试”,即创建一个具有交替正负异常的合成真实模型,从中生成合成数据,然后看反演是否能恢复该模式。虽然视觉上很吸引人,但这可能具有危险的误导性。棋盘格模式,以其规则、与网格对齐的结构,可能恰好与反演的“最佳点”对齐,使得分辨率看起来比对于任意形状异常的实际情况要好得多。它可能会掩盖涂抹是高度各向异性的事实——例如,一个异常在东西方向上的涂抹可能比南北方向严重得多。仔细检查PSF,虽然计算繁琐但科学上严谨,为我们的层析成像视野提供了一份远为真实和详细的“眼镜处方”。

此外,地球的某些部分可能对我们的震源和接收器网络来说是完全不可见的。如果没有射线穿过某个特定区域,该区域就处于问题的​​零空间​​中——数据中完全不包含关于它的任何信息。正演算子的奇异值分解 (SVD) 为我们提供了一种理解这一点的强大方法,它优雅地将模型空间划分为被数据良好约束的部分和完全不受约束的部分。当我们引入正则化——比如我们已经见过的吉洪诺夫阻尼——我们实质上是在为模型中数据无法看到的部分提供一个合理的猜测(例如,结构是“光滑的”)。

超越简单图像:新物理学与更广泛的联系

我们的旅程并未止于标准层析成像。地球比我们最初想象的要复杂,我们的数学工具也比我们想象的更具普适性。

简单层析成像中的一个关键假设是地球是各向同性的——即地震[波的传播速度](@entry_id:189384)与传播方向无关。但这通常是不正确的。地幔岩石中的矿物可能由于地幔流的巨大应变而排列,形成类似于木材纹理的“纹理”。这被称为​​各向异性​​。在各向异性介质中,波速取决于方向。这也导致了​​相速度​​(波前传播的速度)和​​群速度​​(能量传播的速度)之间一个引人入胜的分歧。我们测量的走时受群速度控制。忽略各向异性并使用简单的各向同性模型来解释数据,可能会导致显著的伪影,错误地定位异常或制造假的异常。将各向异性构建到层析成像问题中是现代地球物理学的一个前沿领域,需要与波传播的基础物理学建立更深的联系。

也许,从研究层析成像中得到的至为深刻的教训是,我们认识到其基本思想并非地球物理学所独有。考虑一个来自机器人学领域看似不相关的问题:​​同步定位与建图 (SLAM)​​。一个机器人——比如火星车或自动驾驶汽车——在未知环境中导航。它使用传感器(如激光或摄像头)来构建周围环境的地图,但同时,它必须利用该地图来确定自己的位置。地图帮助机器人定位,而定位又有助于改进地图。

这个问题在数学上表述时,与层析成像惊人地相似。未知的机器人姿态和地图路标是模型参数。它们之间的测量值(里程计、回环闭合)是数据。目标是找到最能拟合所有测量值的模型参数集,这是一个经典的非线性最小二乘问题。其数学结构——一个大型状态向量上的稀疏约束网络——是完全相同的。用于解决它的数值方法也是相同的:我们在层析成像中使用的强大的迭代方案,如高斯-牛顿或列文伯格-马夸尔特算法。甚至确保稳定性和处理病态问题的精妙策略,在这两个领域中都能找到直接的对应物。

这是科学原理统一性的一个优美范例。反演理论的抽象框架为解决从地球核心成像到引导另一颗行星上的机器人等各种问题提供了一个通用蓝图。通过学习层析成像,我们学会了一种思考世界的方式——一种从稀疏和带噪声的数据中推断复杂系统隐藏结构的方式。它是一种工具、一门技艺,也是一种视角,无论我们的科学好奇心将我们引向何方,它都将对我们大有裨益。