地球物理反演

地球物理反演（讲义）主讲：朱良保符号规则：黑体小写拉丁字母代表矢量，黑体大写拉丁字母代表矩阵（非一维）。带下角标的非黑体拉丁字母代表分量。1．前言物理科学中一个非常重要的研究领域就是如何由观测数据来推断物理参数。如：太阳的内部结构，储油层的深度，Moho面的深度，核幔边界的形态等。如果给定物理体系的参数，一般来说，由物理定律能够计算出与观测数据相对比的理论数据。由物理定律根据给定的物理参数计算出数据是正演问题。如图1正演反演问题是根据一组观测数据来重建物理模型。需要强调的是，任何形式的反演过程都必须题，在理想的情况下，运用反演理论可以唯一地求出物理模型。如地震学中的H借助正演手段。没有理论上的正演，就不可能把观测数据有效地与物理参数联系起来，反演就失去方向。对于有些问erglotz-Wiechert反演理论，应用理想的走时数据，假定地层的速度是一维的并随深度单调增加，可以唯一的反演出地层的速度结构。在这个理论中，有两条很重要的假设，其一是数据是无误差的，其二是速度是一维并单调增加的。这两条在实际过程中都不可能达到。尽管非线性反演方案在数学上看起来很漂亮，但在实际运用中却非常有限。首先，准确的反演技术一般只能运用于理想的情况，但在实际中很难满足反演条件，如Herglotz-Wiechert反演理论。另外，准确的反演技术在反演中经常出现不稳定。更重要的一点是，实际物理模型一般是连续的。也就是说，模型的自由度是无穷维的。而实际过程中，观测数据是有限的并且是有误差的。有限的带有误差的数据不足于保证反演的唯一性。BackusandGilbert(1967，反演数据d模型m图1.正反演的传统定义11968)证明了，在线性反演中由于数据量的不足以及误差，反演的不唯一性是必然的。非线性反演更是如此。反演的非唯一性特征意味着，存在许多反演模型能够解释观测数据，由观测数据反演得到的模型不一定就是真实的模型。由图1表征的反演过程过于简单了。我们必须做一些其他的事情。实际的反演过程分两步进行。假如用m表示真实的模型，d表示数据。第一步由数据反演推断出模型m~,这一步为推断。既然有许多模型都能解释拟合观测资料，那么推断出的模型到底与真实模型存在什么样的关系，多大程度上表征了真实的模型，推断出的模型的误差到底有多大,这些问题我们必须回答。没有模型误差分析以及不讨论模型的分辨率，反演出的模型就没有太大的意义。于是乎我们必须作第二步工作，评估模型。所以实际的反演过程可以表征为：反演=推断+评估。图2是实际反演过程的示意图。正问题数据d推断模型m真模型m推断问题~评估问题图2.反演过程是推断加评估一般来说，观测数据是离散的。而模型可以是离散的，也可以是连续的。就模型而言，反演分离散方法和连续方法，处理上有所不同。模型参数与数据的关系有时是线性的，有时是非m),，数据矢量为ddd),,,(L=d，连接数据与模型的矩阵为，称为理论算子。因为它包含了正演计算中所有的物理数学信息。。在实际观测中总会有误差，数据的误差用误差矢量表示。模型矢量与数据的关线性的。对于线性问题，目前的有比较成熟的解决方案。非线性问题比较复杂，还没有找到很好的方案解决模型评估问题。一种方案是对非线性问题做局部近似使其线性化，然后采取循环迭代，逐步接近非线性问题的解，其结果依赖于初始模型的选择。另一种方案是模型空间的全局搜索。目前，无论哪种方案都不能很好地解决非线性反演问题。非线性方法的研究是一个挑战。2．离散线性反演如果模型可用有限个参数来描述，则反演问题称为离散型反演。设模型矢量为Tn21m21作用到模型矢量上为系可以表示为εAmdATmm,,(L=mAεAm+=（1）模型参数的选择带有的一定主观性任意性。例如，描述地球的密度你也许选择地幔是均2匀的度模型就可。如果想更精确一些描述，你可，地核也是均匀的，那么密用两个参数表示以用球谐函数在球面上展开，深度上用多项式来表示，这样就需更多的参数来描述模型。有时对同一个反演问题，不同的人根据自己的喜好用不同的参数来描述模型。由于这种任意性，用模型参数描述的模型就不一定是真实的模型。毕竟真实的模型我们是不知道的。不管用哪种模型参数来描述，总有它的合理性，并对众多模型的选择是一种约束。我们不妨把m叫做理论真模型，或简称真模型。由数据反演出的模型叫估计模型，用m~表示。由于数据的误差及某些我们不清楚的因素，一般来说估计模型不会等于真模型。求解方程（1）的方法很多[如，Menke,1984;Tarantola,1987;Parker,1994]，归根到底是求逆算子，使得g−就=dAm~（2）一般来说，个数是不相同的，所以不是方阵，也就不存在逆矩阵。但广义逆是能够给出的。以后再讨论广义逆怎么求。由（1）+=算子g−A称为矩阵A的广义逆。数据个数与模型的A和（2）可以得出估计模型与真模型之间的关系。把（1）代进（2）得gg−−εAAmAm~（3）矩阵AAg−称为分辨矩阵或分辨核。它是一个算子，表为AARg−≡（4）如果m~IR=，则的分量与的分量一一对应，估计模型完全分辨。否则，估计模型就是（3）式m真模型的线性组合。模型的分辨由数据量和数据的结构决定，而与数据的质量无关。右边的第二项表示数据的误差对估计模型的影响。虽然我们不知道数据误差的细节，否则我们就可以把他们从数据中消除掉，但通过统计分析可以得到数据误差引起的模型误差。设数据的各分量间相互独立，第i个分量的标准误差为idσ。由（3）式得εΔ=Δ−gAm~εΔ=Δd由（1）式知dAmΔ=Δ−g~∴()TgTgT−−ΔΔ=ΔΔAddA估计模型的协方差矩阵为mm)()~(~则第α个模型分量的方差为jgijTigmααασ)()()(2−−ΔΔ=AddA()∑=−=12βαββσdgA（5）存在很大的数大。这是我们不希望看到的。以下讨论广义逆的求法。模型的方差不仅与数据的质量有关还与数据量及由（5）可知，如果广义逆的分量中，则估计模型的误差就可能很数据的结构有关。3讨论：从（4）看，模型的分辨与数据无关吗？3.5最小二乘估计设数据分量的个数大于模型分量的个数。由于数据误差的存在，一般来说线性方程组，严格的表示应该是Amd=是不成立的εAmd＋=。但在实际观测中我们并不知道准确的误差，所以Amd=是超定的矛盾方程组，严格的意义下无解。但是我们可以求最小二乘m意义下最大限度拟合数据的解~。先来看一个简单的例子。由两个物体，其质量分别为1m和2m。称第一个物体的重量是克。称第二个物体的重量是2千克。两个物体一起称得1千重量为2千克。方程为（6）（7）很明显，方程是矛盾的超定系统。因为不可能有第一个物体的重量是1千克，第二个物体的重量2千克，两个物体一起的重量还是2千克。其中一定出现了测量误差。在不知道哪一次测一点，表明方程系统不和谐。反演问题归结为以某种方式来协调这一方程组，求出一个最大限度拟合观测数据的解。⎪⎩⎪⎨⎧==+====2213212211dmmdmdm矩阵A由下式给出⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=110101A量错误的情况下，我们不可能去除某一个测量数据或去除某一个方程而偏重另两个方程。图3是这一问题的图示。方程组（6）在),(21mm平面上对应着3条直线。三条线不可能相交于4图3.线性方程组(6)的几何解释通常的方法是求出一个模型m~mAd~−的范数最小。即2L，使得2~minmAd−（8）满足条件（8）的解叫最小二乘意义下的解。以下来求线性方程组εAmd＋=（9）m~。把（8）写成下角标形式满足条件（8）的模型)~)(~(kikijijimAdmAdS−−=求导数0)~(2~=−−=∂∂kikijlijlmAdAmSδ推得kikTliiTlimAAdA~=设存在，则1)(−AATdAAAmTT1)(~−=（10）5所以方程组（6）的最小二乘解为⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛⎟⎟⎠⎞⎜⎜⎝⎛−−=22112111231~m（11）32~1=m35~2=m图3中的小黑方块就是这一解，他离三条直线的距离最近。由（10）可知方程组（9）的最小二乘意义下的广义逆为TTgAAAA1)(−−=根据（4），分辨矩阵（12）IAAAAR==−TT1)(（12）告诉我们，如果存在，则最小二乘解为完全分辨。1)(−AAT2．2最小范数估计Amd=设模型参数的分量个数大于数据分量的个数，并设方程组是和谐的。在此情况下，不存在。由于方程数小于未知数个数，方程的解有无穷多个。举个例子，两个物体和的总重量为2，求，的重量。系统方程为1)(−AAT1m2m1m2m221=+mm（13）则()11=A⎟⎟⎠⎞⎜⎜⎝⎛=1111AAT其逆不存在。图4是方程（13）的解的图示。很明显在直线上的任意点都是方程（13）的解。6此类问题属于不定方程组解的问题。不定方程组的求解不能用最小二乘准则，因为方程组是图4.方程(13)的解的几何解释和谐的，0mAd＝~−。一般来说，不定方程组有无穷多组解。我们的准则是求最小范数2Lm~m~2minmAmd=满足条件的解，即求方程组的解。求2mAmd=的极小值，满足条件。属于条件极值问题。根据拉格朗日参数法，令)(2AmdΛm−+=TS(14)其中(15)为拉格朗日jijiiimA−)(21nTλλλL=Λ乘数矢量。用下角标表示(2dmS+=λ)求导数并令其为零得ikikkAmmSλ−=∂∂2ikikAmλ=2即（16）进而得ΛAmT=2ΛAAAmT=27（17）ΛAAdT=2设矩阵的逆存在，由（17）得TAA（18）dAAΛ1)(2−=T将（18）代入（16）得dAAAm1)(~−=TT（19）（19）即为模型参数的最小范数的解。将其应用于（13）得2L121==mm（20）矢量的范数范数的概念是绝对值的推广，一个标量，用⋅表示。设ba,是相对于数域F的线性空间中的任意矢量，矢量的范数具有如下性质1．等号才成立仅当,,00aa=≥2．aaαα=，F∈α3．baba+≤+n维实空间nR中的任意矢量a的pL范数定义为ppnLap11⎟⎟⎠⎞⎜⎜⎝⎛=∑=ααa由（19）得广义逆（21）1)(−−=TTgAAAA分辨矩阵为（22）AAAAR1)(−=TT2．3混合问题在运用最小二乘准则时我们假设数据量大于模型参数，方程是超定的矛盾方程组，并且存在。在求最小范数解时，我们假定方程组是和谐的，但没有足够的数据量来1)(−AAT8确定模型参数，方程组是不定方程组，但存在。但是在实际问题中，更一般的情况是，由于误差的存在对于某些模型参数我们有矛盾的数据，而对于另一些模型参数我们又没有足够的数据对其进行估计，并且和都不存在，问题是病态的。即使矩阵的逆存在，方程组也可能是条件病态的，就是说，数据矢量的微小变化引起模型参数的巨大变化，数据的误差在模型空间中被放大了,解是不稳定的。对于这些情况，之前定义的最小二乘准则和最小范数准则都不能解决。我们必须另辟蹊径解决这些问题。Levenberg（1944）提出了阻尼最小二乘法。从数学的观点考虑，病态的结果是矩阵的零或接近零的奇异值引起的。1)(−TAA1)(−AAT1)(−TAAA设矩阵M有本征值nλ和本征向量nvnnnvMvλ=（23）)(IMγ+)(γλ+n则矩阵的本征值是nnnvvIM)()(γλγ+=+（24）这意味着，矩阵的本征值可以通过在其对角线元素上增加一个正的值而加大。这一特性被用来定义阻尼最小二乘解dAIAAmTT1)(~−+=γ（25）由于矩阵的本征值大于等零，如果AATγ的值大于零，则矩阵的零或接近零的本征值的影响就被消除了。)(AAT证明AAT的本征值0≥证：nnnTvAvAλ=Q，则2nnnTTnvAvAvλ=0)(222≥=nnnnTnnvAvvAvAv＝λ，证毕。由（25）表征的解可由求下式的极小值获得22~~mmAdγ+−=S（26）以下脚标形式表示22~)~(ijijimmAdSγ+−=求导数并令其等于零得0~2))(~(2~=+−−=∂∂kikjijikmAmAdmSγ9dAmIAATT=+~)(γ即dAIAAmTT1)(~−+=γ所以（27）