第五章 走时数据的反演.doc

第五章 走时数据的反演在上一章我们根据已知的速度结构对射线追踪和走时曲线的计算问题进行探讨，推导了波的速度只随深度变化的一维速度模型的射线追踪的表达式。三维结构的射线追踪虽然较复杂，但它遵循的原则是类似的。现在我们来研究根据一些观测走时数据，通过反演得到可解释走时数据的速度结构的问题。像人们所设想的那样，反演比前面的问题复杂得多。在全球和地壳结构的研究中，地震学家采取的策略一般是把问题分为两个部分：（1）根据所有可用的数据确定一维的“平均”速度模型。这一般是一个非线性问题，但因为我们要求解的是随深度变化的维一维函数，故比较易于处理，分析没有超出这个范围。（2）如果有足够数量的三维射线的覆盖，把一维模型作为参考模型，对每个数据，把观测的时间减去预测的时间，计算走时残差。对相对于参考模型速度扰动的这些走时残差进行反演就得到三维的模型。如果速度扰动足够小，就可以把问题线性化，即使对大量的数据组，都是可计算的。这是层析反演方法的基础。现在依次考虑上述的每个问题。我们把地震定位问题推迟到这章的最后再讨论，现在假定震源位置已经精确地知道。5.1 一维速度反演在开始反演前有必要设想一下怎样根据走时得到一维的速度结构。假定我们有一个没有三次往返、或没有低速区的简单的走时曲线。曲线上的每一点有一斜率，它给出了射线在转折点处的速度。于是，我们必然知道存在一个特定的速度，问题是这是何处的速度。这相当于确定走时曲线上每一点相应的深度。为此，我们必须知道所论及的深度以上的速度结构。在地表，这是已知的，这使我们的方法可以向下做下去。我们知道在起点（震源）的深度（零）和速度（曲线的斜率）。然后，我们可以研究曲线上与其邻近的点，计算这一点的速度，找出预测的走时曲线通过观测点所相应的深度。按这方法，我们可以沿曲线连续做下去，直到深部。然而，这未必是一个严格的方法，有若干没有回答的问题。总是能得到一个速度模型吗？是否存在不止一个的速度模型，可预测给出同样的走时曲线？现在我们来对这些问题进行讨论，这通常可仿效Aki和Richards（pp:643-651）的论述，读者可更详细地查阅这些论述。回顾一下第四章关于一维速度模型的从地表地表的走时和距离的公式： （5.1）假定有给定的完整的曲线。通过测量曲线的斜率，可以得到，从而得到。我们的目的是反演。上面的反演问题与Abel在1826年解决的很老的问题相似。Abel的问题是查明斜坡的形状，给出球沿斜坡上滚，然后再返回的长度有多大的测量结果，结果是球的初速度的函数（用了早期的物理学问题的不真实的假定：球没有摩擦，没有转动惯量，粘附在斜坡上，不会弹起来）。可以按动能与势能相等，根据初速度计算球可能达到的最高点。Abel阐明了可根据积分变换得到解： （5.2） （5.3）这里是球的最高点，是走时。用作为积分变量，可以把（5.1）式的方程表达为类似的形式： （5.4）这里是时的慢度。现在将其与（5.2）和（5.3）式比较，令，和，即得到： （5.5）分部积分，得到： （5.6）1903-1910年就有三个研究人员在地震学研究中独立推导出方程（5.5）和（5.6）式，称为赫格劳兹维歇尔贝得曼公式（往往称为赫格劳兹维歇尔公式）。对球状地球，可推导出类似的公式（参见Aki和Richards，P.648）。为了用这方程来得到速度的深度函数，我们选取慢度的值。积分上限表示射线参数的射线的射程，是根据曲线得到的。然后对，从0到积分（注意在（5.6）式中是的函数）。这就给出了依赖于慢度的深度。以不同的值重复作这样的计算，就可以得到，从而得到所想要的速度剖面。当有低速区时，这些公式是有问题的。在这种情况下，不连续，没有唯一的解。对这种情况，传统的做法是设想在低速区存在若干有不同速度的均匀层。由于没有射线在这些层里转折，因此低速区的这些层可以任意移来移去，穿过低速区的射线的积分的走时和距离不会改变。用这方法也不能唯一地确定低速区的厚度，而是把它的极限作为最大的厚度。不论解析多么完全，在地震学中很少用赫格劳兹维歇尔（HW）公式，这至少有两个原因：首先，HW假定我们有已知的连续的曲线，而实际上我们总是只有有限的走时点。这意味着走时曲线必须在这些数据之间进行内插，这些内插方案的差距导致于有不同的速度剖面。实际上，会有无数个稍有不同的速度模型，它们都适合于有限数目的点。然而，更严重的问题是实际的地震数据一般多少有噪音，自身相互矛盾。图5.1展示了实际数据的典型的例子。在左边的例子中，小的时间变动导致点的波动，不可能把这些点与能得到期望的1-D速度模型的光滑的、物理上可靠的曲线相联系。在右边的例子中，我们注意到可以把几次地震的数据结合在一起，用这些数据确定“平均”的速度剖面。 图5.1 呈现离散的走时观测结果往往使速度剖面的反演复杂化由于这两个原因，HW公式不能直接应用。下面我们对地震学家反演这些不完善的数据的一些方法进行讨论。HW公式的主要用途是阐明精确的曲线可给出速度剖面的唯一的解。因此，许多反演的策略是把寻找最佳速度模型的问题简化为寻找最佳的曲线的问题。5.2 线性似合速度反演的一个最简单的方法是用一系列直线来拟合走时数据。在50年代和60年代早期，地震学家广泛用其来解释他们的海域折射实验结果。通常是用人工方法，以一定的斜率画直线，其位置是用眼睛估计的（在计算机出现前，这是普遍采用的）。由于直线的斜率确定了地震波的速度，所以可以直接反演数据来确定由少量均匀层组成的简单的模型（图5.2）。通过把每个线段变换为的点，很容易得到这模型。回顾一下多个均匀层情形的公式（4.25）： （5.7）图5.2 用数据似合反演“层饼”速度模型注意，顶层的慢度由第一个线段的斜率确定，第二层的慢度由第二线段的斜率确定，等等。第一层的厚度由前两个线段的斜率和时差确定，即当给出了和，我们就可能解方程得到（也使用，），然后用和解方程得到，等等。通常用三条线段来拟合海洋地壳的走时数据，导出了三个均匀层的模型：层1为沉积层，层2为上地壳（2公里厚），层3为下地壳（4公里厚）。虽然现在我们知道上地壳（层2）包含若干构造，速度不是常数，有很陡的速度梯度，但这些层的名称还是描述海洋地壳的标准术语。对“层饼”模型要注意的是，例如它预测走时曲线有三次往返，在走时曲线的每个拐角附近有第二波至。存在许多可供选择的模型，它们都给出了相同的初至（时间），只是走时曲线的第二个分支的形状有差别（图5.3）。对图5.3底部所示的情况，模型使所有中间的射线聚集在同一距离的地面上（1932年L.B.Slichter首先完成，得到这速度剖面）。如果没有第二支走时曲线的数据，就没办法来区别这些模型。只有我们有一些先验的理由认为有少量的速度近拟为常数的层的速度模型，均匀分层的方法才有意义。图5.3 左边的每个速度模型给出相同的初至（时间），只是走时曲线的第二个分支有差别。5.2.1 拟合曲线的其他方法让我们放弃上述所讨论的简单化的直线拟合方法，直接确定在某种意义下（例如，最小二乘方）逼近所观测到的离散的点的曲线。曲线的任何简单的，离散的描述往往使曲线上有凹形的接缝线，而对第一个初至分支，这在物理上是不可行的。如果我们增加曲线斜率的变化为正值的约束（例如用非负的最小二乘法），就可以避开这一问题。然而，正的约束使曲线有一系列的直线段。起初，这似乎令人惊奇的结果是对反演问题典型的正值约束。曲线在实际上凹的地方相接。从在物理上最接近于数据点的可实现的曲线来说，这是对的最佳拟合。然而像上面所看到的，用一系列线段把曲线参数化难以得到速度与深度关系剖面的唯一的解。而且，即使下面的速度分布是平滑的，分段的曲线也只是一种近似的结果。在我们预期有平滑的速度剖面和没有三次往返的走时曲线的情况下，用平滑的曲线来拟合数据，可以得到最佳的结果。这可以用多项式或样条表达式来实现。走时曲线必须足够弯曲，使得线段上凹不成问题。在三次往返普遍不存在的下地幔，这种近似的应用获得某些成功。然而，由于不同的曲线拟合方法导致不同的曲线（和不同的速度模型），因此不唯一性仍是一个重要问题。因为速度是的非线性函数，因此在反演时，赫格劳兹维歇尔公式使用起来不方便。由于可能有线性公式，在域里做工作可能比较方便（例如Garmay等人，1979）。回顾一下地表地表的射线几何学的公式： （5.8）现在把积分变量从变为，如果单调减小（没有低速区），那么积分限的变换就比较简单： （5.9）这里是在地表的慢度（最大的慢度），在转折点的慢度等于射线参数，如果对这方程作分部积分，即得到： （5.10）不论在上、下限，左边的项都为零，故可写成： （5.11）注意相对于的变化，此表达式是线性的。如果函数乘以2，那么函数也加倍。这种线性关系使线性反演理论的许多方法可用于寻求速度深度剖面的问题。作为一个简单的例子，考虑这样一种情况；已知一系列的射线参数的离散值，我们按一系列的均匀层把模型参数化。的积分方程（5.8）式就变为求和： （5.12）这里是第层的厚度。我们依此写出： 等等，在矩阵形式中，这变为： （5.13）或简写为： （5.14）这里是上述定义的矩阵。注意所有的、和值是已知的，只有包含在矢量里的层的厚度是未知的。因为这是一个线性系统，因此我们可以用标准的方法来求解。如果层的数目少于值的数目，这一般是一个超定问题，可以用最小二乘法来求解。如果层的数目大于值的数目，问题是不确定的。为了得到解，必须有某种形式的规则。如果n=m，我们令，等等。那么就可以对适合于由一系列直线段组成的曲线的层饼模型进行求解。在这种情况下，的第一项和的最后一项被删除。类似于方程（5.7），方形的矩阵有非零的对角线。有理由把这例子进行推广来解比均匀层更平滑的速度模型。例如，可以用与模型每个点（）相联系的线性慢度梯度把速度参数化。在这种情况下，为得到平滑的速度梯度，所需要的模型点的数目要少得多。函数的这种单值的性质（把的三次往返“分开”）和与速度关系方程的线性特性，使在域里做工作比在域里容易些。遗憾的是，地震数据通常是以形式，因此在进行分析前，必须把它变换到域里。原则上，当给定了没有噪音的连续数据，很容易根据计算，反过来也是这样，所用的关系为： （5.15）注意，的一条线对应于的一个点，而的一条线对应于的一个点（图5.4）。图5.4 域里的线相应于域里的点，反之亦然。实际上，我们通常有一系列点的数据。根据一系列的点建造曲线的一个方法是把每个点变换成一条直线，在图解上，这些线交点的包络线就给出所希望的曲线（Bessonova等人，1974，1976）。这种方法不能给出曲线的精确的解，为了根据这些线的包络线给出唯一的曲线仍需要有某种形式的拟合方法。也可以把波型数据用倾斜迭加方法直接变换到域里，线的迭加就给出图像的各个点。如果数据随距离均匀分布，且到时合理相关，此方法的结果很好。它的优点是不需要对地震图作任何挑选。然而，为了反演模型，仍需根据倾斜迭加图像把曲线参数化。有时可以用这种形式的分析来估算曲线的上限和下限。尽管有点主观因素，但可以根据上述图解的包络线或倾斜迭加图像来测量曲线的上限和下限。另一种办法是通过寻找限制曲线的适当斜率的线来得到对特定的值的的上、下限（图5.5）。由于定义下限有某些主观因素，以及在一些情况下，漏检的走时曲线三次往返的第二波至可能落在假定的的上、下限之外，这方法甚至也不能令人满意。如果只可得到初至数据，那么解释三次往返的问题仍然存在，速度反演的唯一性问题仍然受到严重的限制。图5.5 根据走时数据得到的上限和下限可以用来确定的上、下限的位置和速度深度剖面。的两条线展示了最小和最大的深度，在最小和最大深度可能存在特定的速度，仍然满足走时约束。然而，假定我们可以得到曲线的可靠的上限和下限。在这上、下限范围里大量的可能的曲线，相当于有大量的可能的速度模型。我们怎样在这些模型中作出选择呢？我们将讨论两种方法：（1）为得到上、下限的极值反演；（2）对的平滑约束。5.4 线性程序设计和规则化方法或许我们可以采用的最稳健的方法是寻求相应于限度的的限度。通过询问这样的问题来规定的“限度”。对某个给定速度，可能的最小和最大深度是多少？为了得到明确的答案，我们必须再次把低速区排除在外。用像（5.13）式所示的方程把问题公式化，用线性程序设计原理来得到的最大或最小值（例如，Garmany等人，1979；Orcutt，1980；Stark等人，1986）。其结果是所容许的函数的一个狭窄的区域（图5.5）。这些结果有时被错误地解释为表明在这个边界范围里的任一速度模型都是可取的。这是不正确的。在这个边界范围里人们可能给出的多数速度模型都将给出超过限度的曲线。所展示的边界表示的是对某特定的速度所允许的最大或最小的深度。给出某一速度所相应的最大或最小深度的速度剖面并不是真正的界面。这种极值曲线反演方法因其可对速度剖面作严格的限制，因此引起人们的兴趣。然而，实际上这些限制本身往往被证明太宽，以致于意味着的分辨率只能是很粗的（自然妨碍其应用）。于是存在着使用其他方法的诱惑，这些方法似乎可给出比较精细的分辨率，但这未经检验。以上对部分问题作了讨论。简单的最大和最小深度的限制未必能告诉我们关于我们也许感兴趣的曲线的精细结构。线性程序设计方法给出的精细程度可能大于这些限制的方法，但也难以用可以理解的方式标绘出所允许的模型。极值曲线反演的另一个困难是偏离曲线的“点”实际上是对模型作出限定的一些点。大多数的数据或许似乎是相当一致的，但有少数边远的点肯定使和的范围变宽。总是试图去掉这些点，但处理到什么地方算完了呢？减小主观因素的办法是对的误差模型（例如高斯模型）作出假定，然后解包含统计置信限度的最佳拟合的模型（例如，Dorman和Jacobson,1981）。然而，仍然难以对在计算的置信限度里的所容许的模型进行检验。规则化是取代极值曲线限度的方法，为处理反演问题的不唯一性提供了一种可供选择的办法，这是在许许多多的模型中寻求可最大限度地运用一些模型参数的单一的模型。例如，人们可以探索一种与数据相一致的最平滑的模型。这种方法的优点是在最平滑的模型中的任一结构（例如“粗糙度”）都是真实的，而方法本身正试图将其消除。就地震波速度的反演来说，可以用函数的二次导数来度量模型的粗糙度。如果有可用的数据的话，就可以列出问题的线性方程式（Stark和parker，1987）。这类问题的反演结果往往可以用模型的粗糙与数据拟合差之间的平衡曲线（图5.6）来表达。如果我们不顾及对数据的拟合情况，就可能得到很平滑的模型。然而，如果让数据拟合差达到可能的最小值就导致于模型的粗糙度大大地增大。由于总存在着与数据相关的某些误差，因此，试图对数据进行完美的拟合是没有意义的。通常认为最佳的模型位于平衡曲线“拐角”附近的某个位置，在这里数据有较合理的拟合，模型还算平滑（如果数据的统计不确定性是已知的，可以更定量地寻求最佳模量，但这几乎是不可能的）。图5.6 地球物理问题反演的分辨率往往涉及到模型粗糙度与数据拟合差之间的平衡。平滑的模型往往遇到与极值曲线限度模型同样形式的问题，因其看上去过于完美，以致于仿佛人们做了大量的分析也不可能得到。这当然是问题的要害。如果可以用平滑的模型来拟合数据，那么为什么要提出更复杂的结构模型呢？无论如何，这一告诫是适宜的，人们不能开始就想像把最平滑的模型作为“最好”的模型。通常没有先验的理由期望地球是平滑的，最平滑的模型并不是最可能的模型。5.5 归纳：一维速度反演我们之所以对一维速度反演问题作相当详细的讨论，不仅是因其重要性，也是对地震和地球物理反演问题经常遇到的复杂性问题作出说明。尽管作了大量的努力，根据真实的走时数据确定最好的模型，或最好模型的限度的问题仍远未解决。这有几个方面的原因：（1）难以得到三次往返的第二支走时曲线，走时数据存在固有的不唯一性。（2）走时数据通常有噪音，包含随机和系统的误差。系统误差多数来自于在一维速度反演中没有包含在内的横向速度变化。进一步推进一维走时反演方法进步的努力，会带来多少好处是不清楚的。存在着由于统计量分散，不能用所有可用的约束条件来获悉地球结构信息的危险。通常由地震图可得到的信息比走时多得多。例如，地震波振幅对速度梯度是很灵敏的，可以用来诊断三次往返的存在。现阶段的一维模拟技术一般包含模拟整个波形的合成地震图。在这种情况下，走时反演只是用来作为更精细方法的初始模型。往往可以用其他震相，例如来自间断面的反射波或转换波来分辨比用来折射波所得到的，详细得多的精细尺度的结构。人们的注意力正朝解决三维结构方面转变。在解释地球结构时，一维模型的用途有明显的局限。如果数据离散主要来自于横向不均匀性，那么就不能对寻找“最好”的一维模型的问题作出规定，由速度扰动所造成的走时的系统偏移注定了假定数据没有相关误差的任何统计处理都是不成功的。5.6 三维速度反演与用甚至是最好的一维参考地球模型所预测的走时比较，观测的走时通常显示出某些离散。可通过把观测的时间减去预测的时间来计算走时残差，。负的残差是由于波的传播比平均结构的情况快些，因而早到达的结果，而正的残差则是由于慢的速度结构导致波晚到达的结果。在所选择的距离范围里的残差往往可以用残差分布的直方图来表示。残差直方图的展布可以模拟为两部分之和：（1）由震相读数误差引起的随机的离散。（2）由横向不均匀性引起的系统的走时偏差。三维速度反演方法是解横向速度扰动。这些方法与医学上诸如CAT扫描的成像方法类似，现在把这些方法叫做地震层析。然而，值得注意的是地震速度反演比医学问题复杂得多。这有以下几方面的原因：（1）地震射线一般不是直线，而是速度模型自身的函数。（2）地震震源和接收器的分布稀疏，且不均匀。（3）地震定位不精确，往往随速度模型扰动。（4）普遍存在着读数和时间误差。因此，当地震学家们用类似医学上的层析来表达地球结构影像时，可能出错，这是因为“影像”这个术语意指对结构的相当直接的度量，然而实际上，在地震波速度反演中，为了处理上述的困难，需要作许多模型的假设。给出三维速度扰动相对容易些，较富有挑战性的工作是评估统计意义，粗糙度和分辨率。5.6.1 层析问题的建立假定可以得到一维的参考模型，下一步是把三维速度扰动模型参数化，普遍采用两种不同的方法：（1）把模型分成一些速度扰动均匀的块体。（2）就全球模型而言，可以用球谐函数把横向速度扰动参数化，同时用分层的或多项式函数来描述垂向的变化。作为一个例子，现在我们就体波而言，用图解来说明把块体参数化。考虑一个如图所示的分成一些块体的模型的二维几何图形。对每个走时残差都有一条与之相联系的、连结震源和接收点的射线路径。寻找精确的射线路径就构成了两点之间的射线追踪问题，这可能是一项不平常的工作，特别是对必须在整个三维结构里进行射线追踪的层析模型，更是如此。解两点间射线追踪问题的方法包括：（1）在震源处，逐步改变射线的离源角，寻找能与实际接收位置最吻合的射线路径。（2）为了使射线到达所希望的接收点，对射线路径进行扰动。（3）使用节点的有限差分或者图解理论方法。根据费马原理，我们无须得到正确的射线路径，只要相当接近就足够了。因为走时对射线路径的扰动不敏感。一旦确定了射线的几何形状，下一步是求穿过每个块体的走时（虽然原则上来说，这是简单的，但编制计算机程序可能是项很难的工作！）。沿射线路径总的走时扰动为每个块体的走时与块体内的速度扰动的乘积之和。换句话说，走时残差可表达为： （5.16）这里是穿过第个块体的走时，是该块体里的用分数表示的速度扰动（注意没有单位，当速度快1%时，当速度慢1%时，等等）。假定射线路径和是固定的，是根据参考模型进行射线追踪所得到的值。注意如果一维参考模型包含速度梯度，那么各个块体里的速度扰动是常数，但每个块体里的速度可能不是常数。由于速度扰动对射线路径有影响，所以方程（5.16）是一种近似的表达，这种近似只有在的值很小时才是正确的。 如果我们令没有射线通过的块体的走时为零，那么可以把第条射线路径的走时残差表达为： （5.17）这里是模型中块体的总数，因每条射线只穿过模型中很少一部分块体，多数的为零，故几个走时测量有以下矩阵方程： （5.18）这里数字是通过各块体的各个射线走时的实例。我们可以把方程写为： （5.19）这里用常用的符号表示数据矢量，用表示模型矢量，为根据模型预测数据的线性算子。中的数字是每条射线射过每个块体的走时，通常是极其稀疏的，多数元素为零。就所示的情况而言，走时观测的数目大于模型块体的数目。原则上，这是一个超定问题，适合于用标准的方法求解。（5.19）的最小二乘解为： （5.20）由于矩阵总是奇异的，或是病态的情况，因此在层析问题中几乎没有用过这个公式。有些射线路径几乎是相同的，而有些块体没有任何的射线路径通过。对小的矩阵的情况，可以用像奇异值分解（SCD）这样的线性代数方法来减小这些困难。但较普遍的情况是，太大，以致于不可能直接使用矩阵反演的方法。不论是那种情况，实际的情况是反演问题没有唯一的解地下的块体太多，或不同块体之间的扰动互相耦合。 处理病态的最小二乘问题的普遍方法是参考规则化方法的做法，对问题施加附加的约束。规则化方法的一个例子是阻尼最小二乘的解法，按这种解法，可以用下式来代替（5.19）式： （5.21）这里是单位矩阵，是控制阻尼程度的加权参数。问题的最小二乘解是使下面函数达最小值：这里第一项是数据的偏离，第二项是模型的偏差。通过调整参数，我们可以控制偏离与模型偏差之间的平衡。这些约束增加了反演的稳定性没有射线穿过的块体其扰动为零。异常均等地分配在相同的射线路径穿过的各块体之中。但是，阻尼最小二乘解法未必能给出平滑的模型，这是因为由使函数达最小值而给出的是被最小化的模型，不是粗糙度。邻近块体模型的扰动可以很不相同。通常用拉普拉斯算子来度量块体模型的粗糙度，对二维和三维块体情形，这算子可以用差分算子来近似。为了使达最小值，在（5.21）式中，我们用来代替I： （5.22）这里是适用于所有模型块体的拉普拉斯算子的有限差分的近似。目标块体与邻近块体平均值之差给出了每一行的值。例如，在二维模型的情况下，拉普拉斯算子变为：这里是第个模型点的拉普拉斯算子。在这种情况下，最小二乘反演就是使以下函数达最小值：这里控制了偏离与模型粗糙度之间的平衡。与阻尼最小二乘法比较，这类规则化方法使反演的稳定性有不同程度的增大。所得到的模型是平滑的，但偏差未必是最小的。对没有射线路径穿过的块体将用邻近的块体进行内插，而更冒险的是对模型边缘的没射线路径穿过的块体采用外推的做法。阻尼最小二乘和最小粗糙度反演都既有优点，也有缺点。使用哪种规则化方法最佳应随问题的变化而选择。一般来说，人们对给出尺度长度与块体尺度相当的精细尺度结构的阻尼最小二乘法的解表示怀疑，而对最小粗糙度方法的解，只有当所研究的区域只有很少的数据，导致于出现大幅度的异常时，才表示怀疑。至今为止，我们假定所有数据是等权的，但由于走时残差往往不服从高斯分布，存在许多偏离高斯分布的点，因此对数据作均等的加权不总是解层析问题的好思路。不同的方法都提及这一困难。为了消除一些最大的偏离，往往首先取残差窗。在开始反演前，对相似的射线路径的走时残差取平均，得到射线残差和。在迭代过程中，这些异常数据点影响的权重逐步减小，于是看上去像是用了比最小二乘方法范数更强大的误差范数。5.6.2 解层析问题对“小”问题（模型中块体的数目），可以用诸如高斯简化或奇异值分解这样传统的线性代数方法来得到方程（5.21）或（5.22）的精确的解。这时，一个显著的优势就是可以计算分辨率和模型的协方差矩阵。然而。较普遍的情况是太大，以致于难以作这样的实际计算。例如，横向100块，深度方向20块的参数化的三维模型包含了200,000个模型点，显然，我们不可能直接反演一个200,000乘以200,000的矩阵。事实上，我们甚至不能把这样的矩阵装入我们的计算机内存。于是为了解这些问题，我们必须回到为相当稀疏的方程组而设计的迭代方法上来。幸运的是这些方法被证明对于解层析问题是极有效的，已发现能够相当快地收敛到逼近真实解的近似值。这些迭代方法包括ART反投影，SIRT，共轭梯度和LSQR等（参见Nolet，1987，有许多关于这些方法的详细讨论）。虽然领会诸如（5.18）和（5.19）形式的方程有启发性的意义，但实际上我们很少试图建立的矩阵，而是宁可把作为在预测数据的模型方面起作用的线性算子来对待。在计算机计算时，它以子程序的形式出现。由于迭代过程中，一次只能有效地使用的一行，因此有时叫做行作用方法。这些迭代解法的缺点是不可能计算分辨率和模型的协方差矩阵。作为替代，通用的做法是对理论数据组进行试验。通过假定一个特定的速度扰动模型和按实际的数据用相同的射线路径计算走时异常来生成理论数据。然后，对理论数据进行反演，看在多大程度上恢复试验的模型（图5.7）。这种方法的一个例子是脉冲响应试验，在这试验中，单一的局部异常位于感兴趣的区域，看它能在多大程度上被恢复。经常采用的另一个方法是检测板试验，在试验中，对快、慢交错的规则图像的速度模型进行了试验。在这种情况下，方格图像变模糊的程度随模型的位置发生变化，给出了不同区域的相对的分辨率。图5.7 用脉冲响应（上面）和检测板试验（下面）来估算层析模型的分辨率。在每种情况下，按实际数据，用相同的射线路径给出简单速度模型的理论数据组的走时。然后对理论的走时进行反演，看在多大程度上恢复起始的模型。这些试验能否可靠地显示速度反演的精确的分辨率和解的唯一性不总是很清楚的。脉冲响应和检测板试验，因其都假定异常幅度均匀，数据准确，无噪音，因此可能出错。在实际的层析问题中，数据在某种程度上与噪音相混杂，所得到的速度模型包含幅度变化的异常。在这些情况下，往往只能分辨出幅度大的图像。原则上，这些问题有些可以用对数据重新进行随机采样的方法（例如“Jackknife”或“bootstrap”方法）来解决。然而，这要求反演的步骤重复100次以上，妨碍了对计算结果作深入的分析。在层析反演中，有关评估分辨率的问题还没有完全解决，仍是需要继续探索的一个活跃的领域。5.6.3 其他复杂性在前面的讨论中假定震源位置和发震时间已精确地获知。然而就地震而言，这种情况是很少的。由于定位的误差，存在着可能的偏差。以前通常采用一维的参考速度模型进行地震定位。我们期望改为用给定的三维速度模型来定位。确实，在速度异常与地震位置之间存在着耦合的问题。这问题可通过震源和速度的联合反演（JHV），同时求解震源位置和速度结构来解决。在实际工作中，对大量的反演，往往有一个迭代的过程，在这过程中，首先假定一个初始位置，然后导出速度模型。再用新的速度模型对地震重新定位，再导出新的速度模型，等等。只有当各种各样的射线路径可用来确定每个地震的位置（在下一节将对地震定位问题作详细的讨论），地震位置与速度结构之间的耦合才会大大减小。速度反演的另一个模糊的问题与每个地震台站的浅部结构有关。射线一般在每个台站下方，以近垂直的角度向上到达，采样只局限于地壳最顶部很有限的横向区域。同时，通常得不到台站之间浅部结构的信息，因此在大尺度的反演中，各个台站的时间需用台站的修正量作校正，该修正量为到达该台站的所有射线路径的平均残差。就地震定位而言，重要的是根据范围宽的射线路径得到台站修正量，以便使源于深部速度异常的走时差所产生的偏差效应减到最小。当有大量的不同几何形状的射线，且模型的每个单元都有一些彼此夹角范围很宽的射线穿过时，地震层析工作成果是最好的。遗憾的是，情况往往不是这样，这是因为震源和接收器分布不均匀，至少对全球层析问题，又多限于地表。通常，这导致许多块体只有角度范围有限的射线通过。对这种情况，异常只在射线路径的取向上出现（图5.8）。采用规则化或其他数值方法不可能解决这个问题只有补充增加不同角度的射线路径才能使解得到改善。图5.8 当只有角度范围有限的射线时，速度异常的解只限于与射线平行的取向。在某些情况下，存在这样的危险，即三维速度扰动使震源接收器之间的射线路径显著偏离参考模型的射线路径。如果这些射线路径的偏离足够大，那么费马原理就帮不了我们的忙，我们的结果可能出现偏差。所关心的这一问题可以通过对速度模型进行最大限度的三维射线追踪计算来解决，迭代直到得到稳定的解。这就要求做更多的工作，且一般不能做全球层析问题，这是因为对全球，速度扰动只有百分之几。而对局部或区域问题的层析结果有重要价值，这是因为可能发现局部或区域的较大的速度异常、陡的速度梯度或存在的间断面。最后，射线理论近似本身也可能引进了误差，因为理论上假定走时异常只沿实际的射线路径累积。在实际的地震波长范围里，总存在着对邻近的理论射线路径的结构的某种平均。也存在这样的趋势，即射线弯曲或在局部的慢速异常区周围绕射，这可能把偏差带进层析反演中，使可分辨的特征比高速异常的情况低（Nolet和Moser，1993）。在Nolet（1987）和Iyer和Hirahra（1993）的书中可以找到更详细的有关地震层析的方法。5.7 地震定位根据走时数据确定地震的位置是地震学中的一个“古老”的，富有挑战性的问题，一直是地震学研究的一个重要方面。地震通常按发震时间和震源位置来定义。震源是地震的位置，而震中是震源正上方地表上的一个点。在地震定位中，通常把地震作为点源来对待。对于破裂几十到几百公里的大地震，震源不一定是地震的“中心”，宁可把它看成是地震发生时，地震能量开始幅射的那个点。由于破裂的速度小于波速度，可以不管地震最终的尺度和持续时间，只根据初至到时来确定震源。标准目录给出的地震资料，例如震中的初步确定（PDE）依据的是高频体波的走时，不能把这些发震时间和震源位置与长周期反演结果相混肴，后者往往给出地震的矩心时间和位置，表示整个地震的平均发震时间和位置。四个参数描述了发震时间和震源位置，我们把这些参数称为模型，定义模型矢量： （5.23）现在假定我们有在各个地震台的个观测的走时。为了反演这些观测的时间，以得到参数，我们首先必须假定一个参考的地球模型。那么，对每个值我们都可以计算至第个台站的射程和预测到时： （5.24）这里是算子，由其根据给出在每个台站的预测的理论到时。观测的与预测的理论到时之差为： （5.25）这里是第个台站的残差。在某种意义上来说，我们希望得到的是使观测的与预测的到时的残差为最小的。注意，是地球模型和每个台站位置的函数。更重要的是，是模型参数（除发震时间外）的非线性函数。实际上，对一维速度模型，由于可根据已知的参考速度模型的走时表，在适当的范围里对到时进行内插，因此不是特别难以计算的。然而，走时与地震位置参数的非线关系使反演最佳地震模型的工作大大复杂化。这种非线性的影响基至在均匀速度的水平面里二维定位的简单情况下，也是明显的。从座标为的台站到点的走时为： （5.26）这里是速度。在这方程中。不论是与还是都不成线性的定标关系。这结果表明，我们不能用解线性系统的标准方法来得到解。对一组给定的台站走时，不可能用单一步骤的方法来求最佳的地震位置。 在讨论实际的定位策略之前，考虑如果有大功率的计算机可以利用，我们可以做些有益的工作。在这种情况下，我们可在所有可能的地震位置和发震时间范围里作网格搜索，计算每个台站预测的到时。那么，我们可以求得使预测的时间与观测的时间最一致的特定的。怎样规定“最”一致呢？流行的选择是最小二乘法，即求使下式达最小的值： （5.27）这里是台站数目。把平均的平方残差叫做方差。于是我们试图使残差的方差为最小。在叙述模型时你可能听到一个通用的术语是方差缩减（“我仅用两个参数就得到了50%的方差缩减”或“对原始的数据，他们的模型只给出5%的方差缩减”）。这里我们用方差这一术语对残差的展布作不严谨的描述，它与拟合程序中独立参数的数目无关。在统计上更普遍的形式是用来定义方差。这里是自由度的数目（为减去拟合中的独立参数的数目）。对典型的问题，拟合参数的数目比数据的数目少得多，所以和近似相等。最小二乘法常被用来度量拟合差，这是因为在求最小值的问题中，它使方程有简单的解析形式。如果与之间的拟合差是由中非相关的，随机的噪音所造成的，那么最小二乘法将有助于给出正确的答案。然而，在许多情况下，误差不是高斯分布，在这种情况下，最小二乘法对“远离”的数据给的权重太大（残差为2比残差为1对拟合差的贡献大4倍）。作为一种可供选择的替代的办法，我们用残差的和： （5.28）把拟合差的这种度量叫做标准，当“远离”的数据太多时，它与标准（最小二乘法）强大得多。但通常不用标准，这是因为在方程中绝对值的符号产生了很复杂的问题。作为替代像这样的强大标准的一种办法，在最小二乘问题中可以用迭代程序对残差进行加权，在下一步中减小远离的那些点的影响。当然，在对震源位置作大口径的网格搜索时，我们可以直接用所叙述的任一标准。一旦我们规定了拟合差的度量，就可以按最小的拟合差得到最佳的。下一步是估算我们定位结果可能的不确定性。从最小值附近拟合差函数的特征可以看到这些不确定性的某种显示。在我们的二维例子中，如果假定发震时间是已知的（因为是发震时间的线性函数，对给定的某一位置，确定最佳的发震时间是项日常工作），画出作为和的函数的等值线，显然，当远离最小值时，如果迅速增大，我们可得到地震位置比增长缓慢所相应的位置更加精确的解。我们怎样把这个论点定量化呢？因为可以很好地获得高斯过程的统计量，因此最普遍的方法是以最小二乘法和标准为基础，定义： （5.29）这里是期望的由测量误差引起的第个台站残差的标准偏差。期望的的值接近于自由度的值（在我们的情况中，因为有4个元素，所以），通过查阅标准偏差表（例如表5.1），可以得到95%置信限度的值。例如，如果我们用14个走时数据进行定位，那么，在对震源位置作最佳拟合时，根据残差计算的值在3.9418.31之间的概率为90%。值超过18.31的概率只有5%。远离最佳的拟合位置，值将增大，在勾画值的等值线时，我们可以得到相应于这一位置的估算的95%置信限度的误差椭圆。表5.1 分布的百分数（95%）（50%）（5%）51.154.3511.07103.949.3418.312010.8519.3431.415034.7649.3367.5010077.9399.33124.34注意，上述分析中，是临界值完全根据各个走时测量的随机的、非相关的高斯误差所造成的数据偏差的统计量。然而，地震定位问题的拟合差通常比预期的计时和检测的误差大。如果把调整到显著小于平均残差，那么的量度可能表明，应该抛弃这样的解，这很可能是因为对模型速度结果的偏离控制了拟合差。另外，如果被调整到显著大于平均残差，那么因为数据“过于好”，也应该把对震源位置的这种最好的拟合抛弃。为了避开这些麻烦，往往根据在最佳位置的残差来计算数据估算的不确定性： （5.30）这里是最佳拟合位置，对（5.29）式中所有的都用常数值，即： （5.31）注意，使得在最佳拟合震源处的值接近于分布的50%的点。通过绘制的等值线，我们得到解的估算的置信限度为95%的椭圆，也就是说，我们可以逼近真实位置的概率为95%的区域。要注意的是：地震目录通常给出了深度、纬度、经度的标准误差。此标准误差是相应于分布中的第66个百分位数的点。然而，往往把公布的标准误差看成是模型参数各自的扰动，也就是说，它们不一定对误差椭圆的取向可能多少有点斜的事实作出描述。因此，在许多情况下，甚至在统计形式意义上来说，公布的标准误差太小。典型的置信限度的一个更严肃的问题是没有考虑与由于模型没有给出横向不均匀性所引起的走时残差有关的变化。例如，考虑一个垂直断层把地壳分成速度稍有不同的两个块体的模型（图5.9）。由于走时的系统偏差，使得往往把发生在断层上的地震错定为偏离断层，在速度较快的块体这一侧。这种可能性不能用正常的误差分析来解释。在这种情况下，实际上不适当地假定了不同台站之间走时的不确定性是不相关的。假如反演中所用的台站分布很好，那么，没有横向不均匀性模型的影响就是地震定位误差的主要来源。ISC和PDE目录中全球地震定位，其水平位置和深度通常偏离25公里左右（假定用例如pP等深度相位来约束深度。如果没有深度震相，深度误差可能大得多）。改善地震定位的可用方法包括联合速度反演（参见前节）和主事件方法（参见5.7.2节）。图5.9 如果断层两边地震波速度有变化，对发生在断层上的地震，定位往往有偏离。当没有很好分布的台站可利用时，定位误差可能相当大。例如，由于存在着对距离与发震时间之间耦合的问题，发生在台阵外的地震，其距离难以控制得很好（图5.10 ）。在这种情况下，如果有位于地震另一边的一个台站的走时可以利用，地震定位即可得到显著的改善。为了避免定位的这几方面的不确定性，通常要求台站最好围绕地震成最佳方位的分布。另一个问题是当没有近距离的台站可利用时（图5.11），存在着深度与发震时间的耦合，由于射线的离源角很相似，深度的变化可能由发震时间的改变得到补偿。图5.10 对网外地震，定位往往约束得不好。图5.11 如果只有远距离的台站，可能难以确定地震深度在前面的例子中，我们假定只有直达的波可利用。如果在同一台站记录到另外其他的震相，由于使用不同震相的时间差可以消除发震时间效应，因此可使地震定位精度得到显著改善。例如，由传播速度与波不同的波，可直接依据的时间来估算每个台站的震源与接收台之间的距离（按方便的地壳震相的经验法则，以公里为单位的震中距是以秒为单位的时间的8倍）。当用遥测台网数据来确定地震深度时，由于的时间差对深度的变化，反映很灵敏，因此用深度震相可能比用波更好些。5.7.1 迭代定位方法在至今为止的讨论中，我们假定可以直接通过在所有可能的的范围里进行搜索，以求得的最小值