地球物理学中的反演问题.doc

地球物理学中的反演问题1、 介绍物理科学的一个重要的方面是根据数据对物理参数做出推断。通常，物理定律提供了计算给定模型的数据值的方法，这就被称为“正演问题”，见图-1。在反演问题中，我们的目标是根据一组测量值重建物理模型。在理想情况下，存在一个确定的理论规定了这些数据应该怎样转换从而重现该模型。从选择的一些例子来看，这样一个存在的理论假定了（我们）所需要的无限的、无噪声的数据是可以获得的。在一个空间维度中，当所有能量的反射系数已知时，量子力学势能可以被重建Marchenko,1955; Brurridge,1980。这种手法可以推广到三维空间Newton,1989，但是在那样的情形下要求有多余数据组，其中的原因并不是很理解。在一条一维的线上的质量密度可以通过对它的所有本征频率的测量来构建Borg,1946，但是因为这个问题的对称性，因而只有偶数部分的质量密度可以被确定。如果（地下的）地震波速只和深度有关，那么根据地震波的距离，运用阿贝尔变换，这个速度可以通过测定震波的抵达时间来精确构建Herglotz,1907;Wiechert,1907。从数学上看，这个问题和构建三维空间中的球对称量子力学势是相同的Keller et al.,1956。然而，当波速随着深度单调增加时，Herglotz-Wiechert的构建法只能给出唯一解Gerver and Markushevitch,1966。这种情况和量子力学是相似的，在量子力学中，当电势没有局部最小值时，径向对称势只能被唯一建立Sabatier,1973。（量子力学相关概念不熟悉，翻译起来有点坑）图-1尽管精确非线性反演法在数学表达上是美妙的，但它们的适用性是有限的。原因有很多。第一，精确的反演法通常只在理想状态下适用，这在实际中可能无法保持。比如，Herglotz-Wiechert反演假定了地下的波速只依赖于深度并且随着深度单调增加。地震层析成像显示这两点要求在地幔层都不满足Nolet et al.,1994。第二，精确反演方法常常很不稳定。Dorren et al1994已经清楚地展示了Marchenko方程解中这种不稳定性的存在。然而，第三个原因是最根本的。在很多反演问题中，我们要确定的模型是空间变量的一个连续函数。这意味着该模型有无穷多的自由度。然而，在实际实验中，能够用来确定模型的数据数量通常都是有限的。通过变量的简单计算表明这些数据不能承担足够的信息来唯一确定模型。在线性反演问题的背景下，Backus 和 Gilbert1967,1968提出了这一观点，之后Parker1994也提出来这点。这个问题对于非线性反演问题同样相关。在实际实验中有限多的数据可以用来重建具有无穷多自由度的模型这样的事实必然表明反演问题不是唯一的，在这个意义上讲，有很多模型同样可以很好地解释这些数据。因此，从数据反演中得到的模型不一定等于我们想要的真实模型。这意味着图1中展示的反演问题的观点太简单了。对于现实问题，反演实际上包含两步。用表示真实模型，表示数据。由数据我们得到一个估计的模型，这一步称为估计问题（estimation problem），看图2。除了估计一个和数据一致的模型，我们也需要探究估计模型和真实模型具有什么关系。在评价问题中，我们会确定估计模型获得了真实模型的哪些性质以及附带了哪些误差。这部分讨论的实质就是反演=估计+评价。当我们作出一个物理解释却不承认模型中存在误差的事实以及有限的精度，这是没有多少意义的 Trampert, 1998。图-2通常来说，有两个原因可以解释为什么估计模型跟真实模型不同。第一个原因是反演问题的非唯一性，这使得一些（通常是无穷多的）模型满足这些数据。从技术上来讲，这个模型因为模型空间的不充分取样所以零空间存在。第二个原因是实际数据（以及物理理论比我们想要的更频繁）总是受到误差的污染，所以估计模型也受到这些误差的污染。所以模型评价有两个方面，非唯一性和误差传播。模型估计和模型评价对于具有有限自由度的离散模型和具有无穷多自由度的连续模型在根本上是不同的。而且，模型评价的问题只有在线性反演问题上得到很好的解决。因此，离散模型和连续模型的反演是分开处理的。线性反演和非线性反演的情况也是分开处理的。在第2节将讨论有限数量模型参数的线性反演。在第3节中将推广为处理带有无穷多自由度的连续模型的线性反演问题。实际上，很多反演问题都不完全是线性的，但是这些问题常常可以通过做一些适当的近似来线性化。在第4节中将推导出单次散射近似。这种方法形成了运用于反射地震学中的成像工具的基础。Rayleigh原理将在第5节介绍，它是关于线性化的，构成了使用正则模态频率对地球结构进行反演的基础。地震波传播时间层析的线性化方法是基于Fermat原理的，这将在第6节介绍。非线性反演问题要明显难于线性反演问题。第7节将会说明非线性可能是不适定性的一个来源。目前，对于非线性反演问题的评价问题还没有令人满意的理论。在第8节将会介绍三种可用于非线性评价问题的方法。然而，这些方法没有一个是非常令人满意的，表明非线性反演理论是一个有重要研究挑战的领域。2、 解有限的线性方程组在前面的章节中讨论过，反演问题将有限的数据映射到一个模型上。在地球物理学大多数实际应用中，该模型是空间坐标的一个连续函数，因此具有无穷多的自由度。我们暂时忽略这点并假定该模型的特征可以由有限个参数确定。我们将回到这些模型的重要情形，在第3节中这些模型会是无限维的。2.1 线性模型估计对于一个有限维的模型，模型参数可以规定为向量，类似地，数据可以规定为向量。矩阵通过乘积将数据关联到模型上。这个矩阵常常被称为理论算子。确实，在给定的问题上，它包含了我们选择给模型的所有物理和数学信息。实际上，这些数据包含了误差，因此记录的数据和该模型的关系应该是： （1）有一点需要经常注意的是，我们对于包含在模型向量中的模型参数的选择有某种武断性。例如，若想要描述地球的密度，我们可以选择一个模型，在该模型中，地幔和地核具有均匀密度，在这种情况下存在两个模型参数。或者，我们可以把大量定义在球体上的特征方程中的地球密度展开，比如描述横向变化的球谐函数以及描述深度方向变化的多项式，这种情况会有更多的模型参数。在同一个模型上的这两种不同参数化方法对应于不同的模型参数和不同的矩阵。这个例子表明模型m不一定是真实的模型，但是对模型参数的选择通常包含了对于所能构建的模型的等级的限制。以下我们将把认为是真模型，虽然对于它的定义存在很多困难。由记录的数据我们得到模型的一个估计。因为这个估计实际上跟真模型是不同的，我们用来表示估计模型。有很多方法来设计一个逆运算将数据映射到估计模型上e.g. Menke,1984;Tarantola,1987;Parker,1994。无论选择什么估计量，从数据到估计模型之间最一般的线性映射可以写做 （2）算子称为矩阵的广义逆。一般来说，数据的数量不等于模型参数的数量。因此，通常是一个非方阵矩阵，所以它的正常逆矩阵是不存在的。随后我们将说明广义逆矩阵如何来选择，但目前并不需要作详细说明。被估计模型与真模型之间的关系遵循如下表达式（将等式（1）代入等式（2） （3）矩阵称为精度矩阵（resolution kernel），这个算子被定义为 （4）表达式（3）可以写成下列形式来进行解释 （5）在理想情况下，估计模型等于真模型向量：表示我们选择的参数（列在向量中）可以被相互独立估计。等式（5）中最后两项分别解释了估计模型中的模糊度（blurring）和伪差（artifacts）。描述了估计模型向量的元素是真模型向量不同元素的线性组合。我们只能取得模型估计中的参数平均值和模糊度，因为我们无法映射出最完美的细节。在理想情况下，这一项是为零的，此时等于单位矩阵。由（4）可知，对于完美解决的模型参数，精度矩阵为单位矩阵，即 （6）如前所述，通常定义向量的模型参数的定义存在某种歧义。精度算子告诉了我们在估计过程中我们可以独立获得的模型参数的程度。但是，精度矩阵并没有完全告诉我们估计模型和真实的潜在物理模型之间的关系是什么，因为它没有考虑模型参数的选择对于在估计过程中能够得到的模型的限制程度。表达式（5）中的最后一项描述了误差是怎样映射到估计模型上去的。这些误差并不确知，否则它们就能从数据中减去。因为数据中存在误差，所以需要一个统计分析来描述估计模型中的这些误差。当数据不相关且有标准差，则根据数据误差传播，模型估计中的标准差表达为 （7）理想上来看，我们希望同时获得：一个完美的精度，以及不存在误差的估计模型。不幸的是，实际上这是不可能实现的。比如，使用广义逆阵完全抑制了误差传播。这导致（荒谬的）估计模型，这样确实不受误差的影响。但是，这个特殊的广义逆阵对应的精度矩阵是，显然这和理想的精度矩阵相去甚远。因此，实际上我们需要在误差传播和精度限制之间找到一个可接受的平衡点。2.2 最小二乘估计我们现在来考虑这样的情况：独立数据的数量多于未知数的数量。在这种情况下，等式不总是对任意给定的模型都满足，因为数据向量中包含的可能误差使得方程左右矛盾。例如，我们来考虑下面的问题。我们有两个物块质量分别是和。第一个物块的称重得出1千克质量。某人测量第二个物块，结果得出2千克质量。接下来，某人把两个物块放在一起称重，结果发现总质量是2千克。这个问题中测量的结果可以用下列方程组表示 （8）相应的矩阵表示为 （9）显然，这个方程组是不能满足的。不可能第一个物块质量是，第二个物块的质量是，而它们的质量之和。显然测量中存在误差，但是没理由舍弃三个方程中的一个而去支持另外两个。图3（略）生动地阐述了这个问题。在平面中，三个方程对应三条实线。三条线不相交于同一点表示这个线性方程组存在矛盾。所以，采用合理的方法调和这些方程是确定两个物块质量的反演问题的一部分。通常估计模型的一种方法是寻找一个能够最佳拟合数据的模型，在这个意义上，数据向量和再估算数据之间由范数（即欧几里德距离，延伸阅读泛函分析）计量的差值要尽可能小。这表示由模型给出的最小二乘解最小化下列目标函数图-3 （10）详细说明了这个量由下列模型估计最小化 （11）在图3的例子中，最小二乘解是平面中到三条实线距离最短的点，这个点用一个黑色方块来表示。使用矩阵（9），我们很容易得出问题（8）的最小二乘估计量由下式给出 （12）代入数据向量，分别得到估计模型 （13）2.3 最小范数估计在一些问题中，未知量的数量少于参数的数量。例如，考虑这样一个情形：有两个物块和，某人只测量了它们的总质量 （14）相应的矩阵表示为 （15）这个问题被形象地表示在图4（略）中。显然，任何位于该实线上的模型向量都严格满足等式（14）。因此，考虑到质量是正值，则存在无限多完全满足数据的解。一个模型估计可以通过选取一个完全满足数据且具有最小范数的模型来定义，这个模型在图4中用黑色方块表示。图-4对于一个一般的欠定（under-determined 证据不足地说明）方程组系统，最小范数解定义为完全满足数据，即，且最小化的模型。使用拉格朗日乘数法我们能得出最小范数解 （16）给出了详细推导。我们很容易得出系统（14）的最小范数解 （17）2.4 混定问题（mixed determined problems）在最小二乘估计中，我们假定即使由于测量误差导致出现矛盾，我们仍然有足够的信息来求出所有模型参数。所以，这就变成一个完全超定（over-determined 多因素决定）的问题，作为结果是正则的。在最小范数解中，我们假定在可用的信息中不存在矛盾，但是我们没有足够的方程求出所有模型参数。这就是一个完全欠定的问题，这里的是正则的。然而，最一般的情况是我们在一些模型参数上有矛盾的信息，而另一些则因为缺乏信息而无法评估出来。这时，无论还是都不能求逆（invert），这个问题是不适定的（ill-posed）。即使逆矩阵是正式存在的，它们通常也是病态(ill-conditioned)的，这意味着数据向量中很小的变化会导致模型估计中很大的变化。这表示数据中的误差在模型估计中将会被放大。显然，我们需要一点技巧来寻找一个模型，让它对数据中的小变化不那么敏感。带着这个目的，引入一个阻尼最小二乘解。从数学的观点来看，不适定性和病态解是由的零或接近于零的奇异值引起的。（这一句可能不准确）假设我们有一个矩阵，它的特征值（eigenvalue）为，特征向量为，则有 （18）我们容易得到矩阵的特征值为，则有 （19）这表示矩阵的特征值可以通过向原矩阵添加一个乘了系数的单位矩阵来获得。这个性质可以用来定义阻尼最小二乘解 （20）因为矩阵具有正的特征值，则当常数为正时，的特征值将在正方向离零更远（以便消除不适定性和病态）。这样，（20）的解可通过最小化下面的价值方程得到 （21）这个表达式清楚地表明阻尼在当中的作用。最小化（21）式的第一项实际上是寻找最满足数据的模型。最小化（21）式的最后一项相当于寻找有最小范数的模型。一般地，我们不能同时最小化这两项，但是在最小化（21）式时，我们采取折衷的办法去找一个模型既能够合理满足数据又能够获得不太大的模型尺度。参数控制着我们对这两个相互冲突的要求的侧重，因此它被称为权衡系数（trade-off parameter）。对许多应用来说，下面的矩阵性质是非常有用的 （22）在这个表达式中，和是正则方阵，然而并不需要是方阵。这个表达式可以证明当采用加阻尼或者正则化时，最小二乘解和最小范数解（两者都采用了一个阻尼项）是完全相同的。为了理解这点，我们令（22）式中，则有 （23）左手边相当于加阻尼的最小二乘解，右手边是（16）式最小范数解的加阻尼形式。这表明加阻尼后最小二乘解和最小范数解释完全相同的。2.5 最小二乘解的一致性问题最小二乘解看似为寻求超定问题的解提供了一个客观的方法。不过，后面还存在问题。要理解这一点，我们来考虑方程组（8）超定系统。从数学角度看，当我们在最后一个方程左右乘上2，这个方程组系统不会改变。以下这两个方程组系统是完全等价的 （24）原系统矩阵和新的等价系统矩阵可以表示为，以及 （25）在这部分，无撇号的量表示原方程组系统，而有撇号的量则表示变换后的方程组系统。我们容易得到原系统以及变换后系统的最小二乘解（11），如下，以及 （26）代入原数据向量和变换后数据向量的数字值，可以得到下列模型的估计值，以及 （27）问题是，关于同一个模型的这两个估计量是不同的。这很奇怪，因为在（14）式中原方程组系统和变换后的系统在数学上是完全等价的。这两个解不相等是因为原数据空间和变换后数据空间的量度在变换过程中被改变了。这是一种不同的说法，即在使用最小二乘准则解决这两个方程系统的过程中，距离是用不同方法测量出来的。因为最小二乘解最小化了距离，所以当数据空间的量度（或称测量单位）改变后最小二乘解随之发生改变是说得通的。这表明最小二乘解并不像乍看之下那么客观，因为任意变换方程系统会导致不同的最小二乘解！对于最小二乘解，其广义逆表达为。我们很容易推导出原系统和变换后系统的精度矩阵都等于单位矩阵，即：以及。因此，两个系统都具有完美的精度！读者也许有兴趣停下来解释一下（27）式中不相等的估计模型是怎么被调和的。这里出现矛盾的原因是（5）式中的误差传播项。我们知道在表达式（24）中定义的系统中使用的数据一定要带上误差才能使得方程没有矛盾。方程按比例缩放后，两个方程组系统中的数据和误差按不同方式得到调节，这样就得到了不同的模型估计。因此，在方程组中由误差造成的矛盾产生了最小二乘解对任意比例缩放操作的依赖。我们现在来考虑数据向量和模型向量在更一般的变换条件下的最小二乘解。初始方程系统如下 （28）这个表达式不是很正确，因为我们忽略了误差，而误差总是存在的。这就是为什么上面的表达式不能完全满足，我们必须找到这个方程组系统的的最小二乘解。我们来考虑通过一个变换矩阵来对模型参数作变换： （29）以及通过一个变换矩阵来对数据向量作变换： （30）假设可逆，则变换后的方程组系统可以表达为 （31）原方程组系统（28）式的最小二乘解为 （32）变换后方程组系统（31）式的解与上式形式相同，只要令，用替换。这样就给出了变换后模型向量的解。为了将这个解和模型估计（32）式进行比较，我们需要将它变回原模型空间，使用关系式。将最小二乘解写成关于变换矩阵和的形式，有 （33）（表示对先转置后求逆）再假设适当的逆矩阵是存在的，我们对项反复使用矩阵性质，得到，则表达式（33）可以简化为 （34）比较该式和原系统表达式（32）中的最小二乘解，我们发现 当可逆时，对模型向量作变换不影响最小二乘解， 当时，对数据向量作变换也不影响最小二乘解。对于第一个性质的理解，我们回忆前文可知在最小二乘问题中的目标函数不会最小化模型范数，只是对数据做了最小化拟合。对于第二条性质的理解，我们可以对比变换前后被最小的那些量。对于原系统，我们最小化的是： （35）对于变换后的系统，我们最小化的是： （36）当变换矩阵是幺正矩阵（如果一个n阶方阵，它的列向量构成一组标准正交基，那么这个矩阵就是幺正矩阵），也就是当时，上述两个量是相等的。这个结果体现了这样一个性质，即幺正矩阵不影响向量的范数。2.6 最小范数解的调和问题协调性问题不仅针对最小二乘解，最小范数解也面临同样的问题。以（14）中的欠定方程组系统为例，这个问题的最小范数解为 ， （37）作模型向量到新模型向量的一个变换： （38）对于这个新模型向量，方程组系统变成了 （39）需要注意的是，这个变换后的模型向量揭示了这样一个事实，变换后的系统比原系统（14）欠定得更加明显，因为新系统没有对模型参数施加任何约束。变换方程（39）的最小范数解为，。在（38）式变换作用下，变换后模型空间的这个解跟下面的这个原模型空间的解是一致的： ， （40）这个解在图4中用空心方块表示。需要注意的是，这个解与原方程组系统的最小范数解（37）式是不同的。这里出现矛盾的原因跟2.5节中最小二乘解的调和问题是相似的；变换方程（38）已经改变了模型空间的度量标准，所以原模型空间和变换后模型空间中的距离是以不同方式来计量的。一次，原问题和变换后问题的最小范数解是不同我们可以模仿2.5节中的最小二乘解，对模型向量和数据向量在一般变换下得到的最小范数解的变换性质做一个类似的一般化分析。但是，在实际应用中，我们会对方程组作正则化。正如方程（23）所示，阻尼最小二乘解和阻尼最小范数解是完全相同的。因此，一般性的变换性质将在下一节阻尼最小二乘解中来讨论。2.7 寻求更一般的正则化对阻尼最小二乘解变换性质的分析与2.5节中对非阻尼最小二乘解的分析是完全类似的。现在先忽略误差，线性方程组系统按（28）式表示为：，模型向量和数据向量的变换分别按（29）式和（30）式给出：和。再假设存在，则变换后的方程组系统按（31）式表示为：。原系统的阻尼最小二乘解表示为 （41）与（34）式类似，变换后方程组的阻尼最小二乘解表示为： （42）这里的阻尼参数不一定等于原阻尼最小二乘解中的阻尼参数，但是出于我们的目的，我们没必要让它们不等。表达式（42）可以通过（34）式中推导的相同步骤来简化。将项写作，则有 （43）这个表达式指向了一个基本的问题：模型空间的阻尼项由单位矩阵（见（42）式）给出，且当按照原模型向量来表达时，阻尼项为（见（43）式）。这表示阻尼对于模型参数变换并不是不变的。当进行数据向量变换时就会产生项。这表示在对数据向量或者模型向量作变换之后阻尼最小二乘解一般是会变得。因此，需要有一个更一般的正则化可以允许以一致的方式改变模型和数据空间。这样一个一般的正则化可以在（43）式中找到，令，同时定义。这样，一般性最小二乘解表示为 （44）该解最小化下面的目标函数： （45）这个表达式表示一般情况下权阵和可以是任意的（只要它们是正定的以保证有极小值）。按照这种写法，可以被看作是一个权衡参数，它折中了模型的两个特征：模型大小以及模型对数据的不符合性。模型的两个独立的性质不能都取得任意小，因此需要取一个平衡。然而，选取一个最适宜的并不那么容易。我们已经清楚地表明，当你从一个简单阻尼最小二乘解开始时，你可以把那个问题转换成另一个坐标系中的一个更一般的正则最小二乘解，反之亦然。这表明相对于（44）式这种更一般的最小二乘解我们偏向阻尼最小二乘解是没有道理的。实际上，大多数反演问题都是不适定的（部分欠定部分超定）和病态的（数据中的微小误差引起很大的模型变动），它们和零空间同时出现，因此总是有非唯一解。所以需要进行正则化，但是在正则化的选取上存在很大不确定性Scales and Snieder, 1997。这反映出我们在解决反演问题时面临的基本难点：相对于选取正则化，解决方程组系统只是一个次要问题。采用贝叶斯统计法是其中一种方法，我们从统计学的角度结合关于数据和带实测数据的模型的先验信息来讨论反演问题Tarantola and valette, 1982a; Tarantola and Valette, 1982b。权阵反映了我们所掌握的关于数据和模型的真实的物理先验信息（从统计学角度），且独立于测量数据。比如，其中包括获取数据用的仪器的统计噪声特性，也包括来自其他自变量的模型和数据信息。（比如，地球的质量密度一定是正的。）在贝叶斯方法中，权阵为 ， （46）这里的和分布是数据和模型的先验协方差矩阵： （47） （48）在这些表达式中，尖括号代表期望值。在这个解释中，当误差服从高斯分布时，估计量（44）与最可能的后验模型是一致的。如果我们承认先验信息包含真实物理意义这样一个规律，则贝叶斯变换的统计学基础会推导出一个客观解。然而，实际上我们应该意识到对于数据和模型的先验分布的选择常常也是很主观的。读者可以在Scales and Snieder中的专栏“贝叶斯还是非贝叶斯”找到更加详细的信息。另一种方法是定义一个非拟合函数，这个方法赋予模型给定的性质（比如小，光滑等等）Parker,1994。选择先验信息相当于定义一个适当的范数，在这个范数中非拟合的数据和任何给定的模型性质都会被测量。在我们的这个情况下，权阵会按范数定义一个特定的矩阵。选择权阵的一个很常见的例子就是奥卡姆反演Constable et al, 1987，该反演寻求一个与数据保持一致的最为平滑的模型。在这里，我们要求这个模型的梯度平方尽可能的小，也就是（45）式中最后一项是一个离散化，而不是在模型范数上施加一个约束，因此跟拉普拉斯算子的一个离散形式是一致的。2.8 权阵的变化规则反演解的一个基本的要求应该是反演问题的解不依赖于对模型向量或数据向量所做的任意缩放变换。遗憾的是，这个要求经常被忽视，使得不同模型的对比变得毫无意义。了解实际含义可参阅Tramper and leveque, 1990以及Trampert et al, 1992。在这里，我们将推导在使最小二乘解保持不变的变换下，权阵和应该怎样调节。让我们首先考虑（29）式中的模型向量缩放：。在这个变换下，（45）最小二乘量中的模型项将作如下变换 （49）其中 （50）对模型权作这样的变换不会改变最小二乘准则，因此，当模型权阵经过变换后最小二乘解是不变的。值得一提的是，这个变换规则意味着，对于贝叶斯反演（其中的权阵是模型协方差阵的逆矩阵，即），协方差阵应该作如下变化 （51）我们容易从定义式（48）核实这确实是协方差算子的变换规则。接下来，我们考虑数据向量的变换（30）式是怎样影响数据权阵的变换的。矩阵在变换的作用下作缩放变换。在这个变换作用下，最小二乘量（45）中的数据项将作如下变换 （52） 其中 （53）对于贝叶斯反演，数据权阵是数据协方差阵的逆阵（），这样对于一个贝叶斯反演，应该作如下变换 （54）再次注意，这是在（47）式中定义的正确的协方差阵变换规则。这表明贝叶斯观点（即和是模型和数据协方差阵的逆矩阵）确保了在模型向量以及/或者数据向量作变换后，解是不变的。尽管我们已经推导出在模型和数据向量变换下权阵和是怎样变换的，这绝不意味着这些权可以被明确地定义。一个不适定的以及/或者一个病态的反演问题只有当我们愿意通过施加一个正则化（regularization）项来控制解时才能被解决。一般来说，选取权阵和的方法并不唯一。选取这些矩阵则是用户的主观介入。2.9 解决线性方程组系统应该注意的是，最小二乘解总是要求解决一组线性代数方程。比如，方程（44）可以写成 （55）这个公式表示的是线性方程组的一个正方系统，也就是所谓的法方程组。如果我们仅仅是对问题的估计部分感兴趣，那么不需要转换。如果我们还对问题的评价部分感兴趣（事实上我们总是会这样），那么就必须意识到需要被转换，这就要付出额外的计算机运算时间。许多标准的子程序包是可用的，Press et al, 1989对这个问题给出了一个很好很实际的介绍。但是读者需要明白，对于实际的地球物理问题，系统可能变得非常庞大，所以值得去考虑一个专门的程序，这个程序最匹配的性质（对称性，带状的，稀疏的，等等）。法方程组的维数也值得考虑。记住矩阵的维数是，其中是数据数量，是模型参数数量。系统（55）具有模型空间的维数，但是运用（22）式我们可以得到一个与数据空间维数严格等价的系统。选择最小的维数写出法方程可以节省很多计算机计算时间。大多数解决代数方程组的技术直接和作用，总的来说需要足够的存储空间来承载这个矩阵。但是，比如在全球走时断层扫描中，这些维数可能会极度庞大（，），这样就需要采用迭代方法，这种方法一次只能作用于的一部分。另一个经常出现的问题是即便中包含了正则化，它仍然是奇异的或者数字上非常接近于奇异。奇异值分解（SVD）是一个强有力的技术，它能精确判断出问题是什么，并且给出一个有用的数值解。SVD是反演理论中解释获得某个结果最有效的工具。迭代法或者SVD不需要作用于方正系统，因此可能直接使用矩阵。在这部分内容中，明白一般化最小二乘解（44）等价于简单最小二乘解是有用的。 （56）读者可以参阅Tarantola1987中对于正定矩阵平方根的意义的讨论。还要记住在选择权阵时的某种自由（见2.7），读者也许想要直接定义而不是。表达式（56）表明，正则化会往线性方程组系统中添加额外的行，但是扩大的系统仍然保持形式，其中矩阵和数据向量被导致正则化的额外的行扩大了。此时在法方程中，不再是方阵。我们现在将更详细地阐述应用在系统中的奇异值分解和迭代法的本质。2.9.1 奇异值分解（SVD）其奇异值分解是由Lanczos 1961年提出来的，这个方法是将矩阵的特征向量分解一般化至非方阵的情况。我们首先来考虑一个真实的对称矩阵方阵，它的特征向量是，对应的特征值是。对于这样矩阵，其特征向量构成了一个单位正交组，因此任何向量都可以由这些特征向量表示：。当作用于这个表达式时，结果可以写成 （57）用同样的特征向量分解向量，则有，将其代入（57）式，可得下列解向量的展开式： （58）可以看出，小的特征向量可以导致解的不稳定。奇异值分解将这个展开式推广至非方阵矩阵。Lanczos 1961以及Aki and Richards 1980给出了该方法的详细内容。现在来考虑下面的非方阵方程组系统： （59）奇异值分解基于在单位正交组特征向量下的展开式和在单位正交组特征向量下的展开式。这些向量不可能是的特征向量，因为不是方阵，因此它没有任何特征向量。不过，这些向量与下列关系式存在关联： ， （60）容易看出，当向量是的特征向量时，向量是的特征向量，因此这些向量可以很容易被确定下来。和有公共的非零特征值。被称为的奇异值。当作用在上时，结果可以写成： （61）和的上限由非零特征值的数量来确定，因为零特征值对求和没有贡献。数字可以明显小于问题的维数：且。将向量和向量作为矩阵和的列是很方便的，从指标往后（即从开始）的特征项对应于零特征值，需要把它们补充进矩阵和中使得矩阵完整：， （62） ， （63）特征向量的正交性意味着，。特征向量的完整性意味着，。因为特征向量的正交性也存在于和的子空间和中，所以有，。但是，这些子空间中的向量一般不构成完备组，所以一般，。（61） 式一般化至非方正系统可以写成， （64）其中为 （65）从（61）式可知，当作用于一个向量时，对应于零特征值的特征向量不参与贡献。这些特征向量被安排在子矩阵中。这相当于是说根据表达式（64）可知矩阵可仅由和建立。和是没有被算子照亮的空间黑点（我认为这是比喻手法）。因为，所以预测数据与子空间正交，见图-5。图-5这表示中数据向量的任何分量都不能被任何模型所解释。这些数据矩阵的分量必然和数据中的误差或者算子中的误差是对应的，它们是对物理问题的一个描述。因此，被称为数据无效空间。在最小二乘反演中，我们的目的是最小化非拟合数据。最小化非拟合数据实际上是找到一个模型，这个模型产生子空间中的一个最接近真实数据的数据向量。从图-5可知，这只要简单地把包含在数据中的的分量从问题中投射出去。这实际上就是将（64）式中对特征值的求和限定为只对非零特征值求和。当然，当为空时，我们总可以找到那个完全解释数据的，因为是跨越了整个数据空间的。类似地，对特征值的求和限定在非零特征值之上，其产生的影响是模型估计只产生于子空间中，估计模型在中没有分量。在中的任何模型分量都不会影像数据，因为。这表示定义了模型无效空间（model-null-space）。数据跟落在中的模型向量的分量是没有关系的。将模型无效空间中的模型向量分量设置为零意味着在模型估计中我们只能考虑非零特征值。将在向量中展开，同时将在向量中展开，并只考虑非零特征值，我们可以因此用以下的方式将（58）式中的解一般化至非方阵系统中：， （66）使用矩阵和，上述结果也可以写成：， （67）其中， （68）类似于前面的问题，反演问题不是和的函数。如果这两个子空间都为零，算子存在精确的逆。如果存在，我们可以得出残差与是正交的，因此残差是最小二乘准则下的最小值。如果存在，则解（67）在中没有分量，因此是满足最小范数准则的。很明显，当乘以，且奇异值比较小的时候，中的小误差会导致中出现较大误差。这个误差放大过程可以通过限定（66）式中仅对与零相差很大的特征值求和来加以控制（略去接近于零的特征值）。亦或者，我们可以将替换为，其中是一个正的常数。我们可以得出，这与（20）式中的阻尼最小二乘解是等价的。参考Matsuura and Hirata 1982中对这些不同策略的讨论。需要注意的是，去除或者抑制小特征值会导致不同的结果。这使得定量比较同一个问题通过这些根本上不相同的策略所获得的解是几乎不可能。奇