1、第五章第五章 反问题基本方法反问题基本方法内容提要第一节 引言 第二节 地球物理反问题的广义线性逆 有关地球内部，特别是深部的知识，主要还是来自地球物理的观测资料的解释取得的。因此，地球物理学的基本问题就是用地面或地表附近的各种观测资料定性或定量地对地下的地质结构和矿产资源做出判断，这就是解地球物理反问题。用公式表示为 gmGd，在线性常定条件下表示为1mA d。d数据向量，m地球模型的参数向量，gG正演问题dG m中算子G的广义逆，1A为线性算子或矩阵A的逆。第一节 引言 地球物理反问题有广义和狭义两种提法。广义地球物理反问题是指根据地面观测数据不管用什么样的方法和计算手段，能对地下的地质情
2、况做出定量或半定量的估算都可称为解反问题或反演问题。例如合成声波测井曲线的计算、神经网络的地层横向预测、油气储层的多参数分析、矿体埋深与大小的估算和求波动方程系数以及地球物理层析成像、地球物理资料的地质解释等均属于广义反演范畴。狭义地球物理反问题是指与各种地球物理场方程，即各种偏微分方程有关的地球物理反问题的求解方法和计算问题。狭义的地球物理反问题与数学物理反问题或数理方程反问题是密切相关的。第一节 引言 狭义地球物理反问题又可分为两类，第一类是根据现在的物理状态确定物理过程的过去状态。例如用已知物体的当前温度去确定初始的温度分布；根据空间中波场的当前分布计算以前某时刻波场的空间分布等。前者用
3、于地热场的研究，后者用于地震偏移成像、电磁波成像和地震层析成像的研究。第二类是从微分方程的解的某种泛函来求方程的系数或右端项。第一节 引言 第一类反问题求解相对容易，虽然场函数与介质参数二者之间存在不确定性，但计算总是稳定的。因此，在实用中取得了良好的效果。无论是地震偏移问题，还是求地热场的原始分布问题，都是通过解微分方程计算波场或热流场的传播逆过程，虽然这个过程也是物理不可实现的，但是它的物理意义是非常确定的，而且这个过程是唯一的。因此，反演的结果在参数未确知或计算方法非最佳的情况下只存在误差，不存在非解的问题。第一节 引言 第二类数理方程反问题一般是不适定的和非线性的。不适定性表示这类反问
4、题的解不唯一。非线性说明求解上的困难。这类反问题经过半个世纪的求索，至今仍然不是那么好解决的，特别是高维波动方程反问题的研究，无论从理论上，还是从实用上来看，基本没有一个令人满意的求解计算方法。虽然在各个学科领域发表了大量的论文和著作，但一触及实际问题都显得无能为力。特别是计算地球物理反问题由于数据不足，使本来的不适定问题的求解更加困难。第一节 引言 地球模型和数据之间的关系由 dGm （1）给出。假定（1）是线性的。求解线性反问题的最简单的方法，建立在估计的模型参数estm和预测的数据preestdGm的长度或大小的度量的基础之上的长度法。反演的目标是，找出让观测数据与预测数据之间的某种长度
5、达到最小的那一组模型参数estm。目标函数E表示为 11lNlobspreiiliEdde （2）第二节 地球物理反问题的广义线性逆 le表示向量e的 l 范数，1、2 范和无穷范的定义分别为 111:NobspreiiiLdde，（3）1 22221:NobspreiiiLdde ，（4）:maxobspreiiiLdde 。（5）用 2 范作为目标函数，则得到的解成为最小二乘解。此时 111NMMTTiijjiijjijjEdG mdG me edGmdGm（6）第二节 地球物理反问题的广义线性逆 目标函数E对模型参数求导并令其为零，得到 1110,1,2,.,MNNkiqikiqikii
6、mG GG dqM ，（7）将其写为矩阵形式为 TTG GmG d0 。（8）如果1TG G存在，则（8）有解 1estTTmG GG d。（9）第二节 地球物理反问题的广义线性逆 最小二乘解存在与否与dGm所提供的信息量有关。dGm提供的信息不足，不能唯一地确定最小二乘解，此时的反问题称为欠定的。dGm包含的信息太多以致于使其不具有精确解时，则称为超定的。dGm恰好具有足够的信息来确定模型参数，则称为常定的。S R S R S R S R S R S R S R S R S R 欠定 常定 超定 第二节 地球物理反问题的广义线性逆 对于纯欠定问题，未知数的个数多于数据的个数，也即多于方程的个
7、数，此时反问题有无穷个解。为了获得反问题的解estm，必须拥有能精确地选出其预测误差 E 为零的无穷多个解中某一个解的方法。要做到这一点，就得把某些未包含在方程dGm中的信息附加到该问题中。这些附加信息称为先验信息。先验信息的例子：在直线拟合中，它必须通过一个已知的确定点；密度参数必须大于零；地震波速度参数介于 1000m/s-4000m/s；等等。第二节 地球物理反问题的广义线性逆 还有一种先验信息是，反问题的解是简单的，“简单”用解的欧几里德长度度量。此时，反演问题的目标变为：求在edGm0约束下使解的长度TL m m取极小的estm。利用拉格朗日乘子法构造该反问题的目标函数 21111N
8、MNMi iiiiijjiiijmLemdG m （10）（10）式对模型参数分量求导并令其为零，得 2TmG 。（11）第二节 地球物理反问题的广义线性逆 2TmG 。（11）与约束方程dGm 一起求解。把2TmG 代入dGm，有 12TdGG。（12）如果TGG 的逆存在，则可求得拉格朗日乘子 12TGGd。（13）把12TGGd代入2TmG 得到模型向量的估计 1estTTmGGGd。（14）第二节 地球物理反问题的广义线性逆 上述求解线性反问题是建立在预测误差和解的简单性检验的基础上的，得到形如estmMd的解，其中M为一矩阵。反问题的求解实质是M的求解，知道M后就知道了反问题的解。M
9、称为广义逆，用符号gG表示。因此，Gmd的反问题解为estgmGd。最小二乘问题的广义逆为1gTTGG GG；最小长度欠定问题的广义逆为1gTTGGGG。第二节 地球物理反问题的广义线性逆 在某些方面，广义逆类似于普通矩阵的逆，方阵方程Axy的解为1xA y，而反问题Gmd的解为estgmG d。不过这种相似性是非常有限的。广义逆并不是在通常意义下的矩阵的逆。它不是方阵，并且既不要求gG G等于单位矩阵，也不要求gGG等于单位矩阵。第二节 地球物理反问题的广义线性逆 数据分辨矩阵 假定已经找到了一个在某种意义下能够求解反问题Gmd的广义逆，这样就得到模型参数的一个估计estgmGd。这一模型参
10、数估计值对数据的拟合程度如何？把该估计值代入方程Gmd，则得 （15）NN阶方阵gNGG称为数据分辨矩阵。它描述了预测值与数据的拟合程度。如果NI，则preobsdd，因而预测误差等于零。反之，如果NI，那么预测误差就不等于零。第二节 地球物理反问题的广义线性逆 数据分辨矩阵 如果数据向量d的元素具有自然顺序关系，那么数据分辨矩阵就具有一个简单的解释。例如，考虑用直线拟合数据点（z,d）的问题，其中数据是按辅助变 z 的值排列的。如果N不是一个单位矩阵，但是接近于一个单位矩阵（在其最大元素靠近主对角线的意义下），那么从该矩阵的构形就可以看出，能够预测出相邻数据的平均值，却不能预测出单个数据。第
11、二节 地球物理反问题的广义线性逆 数据分辨矩阵 考虑N的第 i 行，如果该行中除了第 i 列上的元素不为零而其余全为零，那么就有可能准确地预测出id。反之，假定该行的元素为 （16）式中的 0.8 位于第 i 列上。则第 i 个数据为 （17）该预测值是三个相邻的观测数据的加权平均值。如果真实数据随辅助变量缓慢地变化，那么该平均值可能产生一个接近于观测值的合理估计值。第二节 地球物理反问题的广义线性逆 数据分辨矩阵 N 的每一行描述了相邻数据能被独立地预测或分辨的难易程度。如果数据具有自然顺序关系，则N 的行元素随列下标变化的图像反映了分辨率的清晰度。如果这些图像具有一个尖锐的极大值，并且其中
12、心在主对角线上，那么数据就得到很好地分辨。如果这些图像非常宽阔，那么数据就得不到很好地分辨。第二节 地球物理反问题的广义线性逆 数据分辨矩阵 因为数据分辨矩阵的对角线元素表示数据在其自身的预测中具有多大的权，所以经常把这些对角线元素挑选出来并称之为数据的重要性n （18）数据分辨矩阵不是数据的函数，而仅仅是数据核G（它体现了模型及实验的几何特征）以及对问题所施加的任何先验信息的函数。这样不用真正进行实验就能计算和研究数据分辨率矩阵，因而数据分辨矩阵可以作为实验设计的重要工具。第二节 地球物理反问题的广义线性逆 模型分辨矩阵 数据分辨矩阵表征了数据是否能被独立地预测或分辨。对模型参数也能够提出相
13、同的问题。想象存在一个真实的、但未知的满足trueobsGmd的模型参数向量truem。该模型参数的某一特定估计值estm与这一真实解的逼近程度如何？将trueobsGmd代入估计的模型参数表达式estgobsmGd中，得到 （19）第二节 地球物理反问题的广义线性逆 模型分辨矩阵 （19）其中的R 是MM阶的模型分辨矩阵。如果R 不是一个单位矩阵，则模型参数的估计值就是真实模型参数的加权平均值。假如模型参数具有自然顺序关系，则分辨矩阵的每一行的图像可用来确定模型中多大尺度的特征确实能够被分辨出来。与数据分辨矩阵一样，模型分辨率也只是数据核和对问题所附加的先验信息的函数，它与数据的真实值无关，
14、因而它可以作为实验设计中的又一个重要工具。第二节 地球物理反问题的广义线性逆 单位协方差矩阵 对于线性函数 mMdv，均值和协方差之间的关系为 mM dv covcovTmMd M 广义逆的协方差为 covcovggTmGd G 第二节 地球物理反问题的广义线性逆 单位协方差矩阵 模型参数的协方差covcovggTmGd G取决于数据的协方差以及由数据误差映成模型参数误差的方式。其映射只是数据核和其广义逆的函数，而与数据本身无关。因而为了描述映射中误差的放大程度定义一个单位协方差矩阵是非常有用的。如果假定数据是不相关的，并且所有的数据具有相同的方差2，那么单位协方差距阵为 2covcovggT
15、ggTumGd GG G （20）第二节 地球物理反问题的广义线性逆 单位协方差矩阵 甚至当数据是相关的，也常常能找到数据协方差矩阵的某种归一化，这样就能定义一个与模型协方差矩阵有关的单位数据协方差矩阵covud，其关系式为 covcovggTuumGd G （21）因为单位协方差矩阵也像数据和模型分辨矩阵一样，与数据的真值和方差无关，所以它也是实验设计中的有用工具。第二节 地球物理反问题的广义线性逆 单位协方差矩阵 作为例子，考虑用直线拟合数据(z,d)。截距1m和斜率2m的单位协方差距阵为 2221coviiuiiizzzNNzzm （22）仅当数据的中心位于z0 时，其截距和斜率的估计值
16、才是不相关的。整个方差的大小是由分数的分母决定的。如果所有的z近于相等，那么该分数的分母就比较小，但是其截距和斜率的方差却比较大（图 a）。反之，如果z的分散性较大，那么该分数的分母就比较大而方差却比较小（图b）。第二节 地球物理反问题的广义线性逆 (a)(b)某些广义逆的数据、模型分辨率和协方差（covudI）最小二乘法 最小长度法 最小二乘解和最小长度解之间存在许多对称性。最小二乘求解完全超定的问题，并且具有完好的模型分辨率，最小长度法求解完全亚定问题，具有完好的数据分辨率。将会看到，求解介于上述两者之间的混定问题的广义逆的数据和模型分辨率介于这两种极端之间。第二节 地球物理反问题的广义线
17、性逆 模型分辨矩阵R数据分辨矩阵N分辨率和协方差的优度评价 与能通过度量模型参数的总预测误差和简单性来定量表示模型参数的好坏一样，可给出几种方法来定量表示数据和模型参数的分辨矩阵及单位协方差矩阵的优度。当分辨矩阵是单位矩阵时其分辨率最高，所以可以根据非对角元素的大小或展布对分辨率进行评价。分辨率展布的优度的度量建立在分辨矩阵与单位矩阵之差的L2范数基础上。有时把这些展布称为狄利赫莱展布函数。因为模型参数的单位标准差是由数据到模型参数的映射中误差放大量的一种度量，可以用它估计单位协方差的大小，即 协方差大小的这种度量未考虑单位协方差矩阵中非对角元素的大小。第二节 地球物理反问题的广义线性逆 具有
18、良好分辨率和协方差的广义逆 考虑是否有可能把这些度量准则作为指导原则来推导广义逆。首先是定义解的预测误差和简单性的度量，然后利用这些度量推导模型参数的最小二乘和最小长度估计。超定情形 考虑形如Gmd的纯超定问题。假设解为estgmGd，并试图通过使上述优度度量的某一组合达到极小来确定gG。因为前面已经注意到超定最小二乘解具有极好的模型分辨率，所以这里试图通过只对数据分辨率的展布取极小来确定gG。现检查N的第 k 行的分布kJ：因为kJ 都是正的，所以使每个kJ 达到极小就能使总展布spreadkJR达到极小。第二节 地球物理反问题的广义线性逆 具有良好分辨率和协方差的广义逆 超定情形 因而把数
19、据分辨矩阵的定义式gNGG代入kJ 的表达式，然后求kJ 对广义逆矩阵的每个元素的极小 。分别对kJ 的三项中的每一项求导数，第一项为 第二项为 第三项为零，因为第三项不是广义逆的函数。第二节 地球物理反问题的广义线性逆 具有良好分辨率和协方差的广义逆 超定情形 用矩阵形式完整地写出该方程便得到 因 为TG G是 方 阵，所 以 左 乘 它 的 逆 阵 而 求 得 广 义 逆1gTTGG GG，该式具有与最小二乘广义逆完全相同的形式。于是既可以把最小二乘广义逆解释成使预测误差的 L2 范数极小的广义逆矩阵，又可以解释成使数据分辨率的展布极小的广义逆矩阵。第二节 地球物理反问题的广义线性逆 具有
20、良好分辨率和协方差的广义逆 亚定情形 在一个纯亚定问题中，数据能够得到精确地满足。因而其数据分辨矩阵刚好是一个单位矩阵，这样其展布就等于零。于是就有可能通过求模型分辨矩阵的展布对广义逆矩阵的每个元素的极小来导出该问题的广义逆。由该方程所 得 到 的 广 义 逆 刚 好 就 是 最 小 长 度 广 义 逆1gTTGGGG。因而既可以把最小长度广义逆解释成使最小长度解的 L2 范数极小的广义逆矩阵，又可以解释成使模型分辨率的展布为极小的广义逆矩阵。这就是最小二乘解和最小长度解之间的对称关系的另一个方面。第二节 地球物理反问题的广义线性逆 具有良好分辨率和协方差的广义逆 包含狄利赫莱展布函数的一般情
21、形 有时要寻找使分辨率的展布与协方差大小的加权和为极小的广义逆，即求 的极小。式中是任意的加权因子。其结果是一个广义逆的方程：从矩阵的代数函数角度来说，从该方程得不到gG的显式解。然而，对于大量特殊选择的加权因子，却能够写出其显式解。例如，当1231,0时，就回到最小二乘解，而当1230,1,0时，就回到最小长度解。更为有趣的情形是1231,0,等于某一个常数（比方说2）及covudI，这时广义逆 此式刚好就是阻尼最小二乘广义逆。阻尼最小二乘广义逆也可以解释成使数据的展布与协方差大小的加权组合为极小的广义逆矩阵。第二节 地球物理反问题的广义线性逆 旁瓣与 Backus-Gilbert 展布函数
22、 如果数据或者模型参数具有自然顺序，那么狄利赫莱展布函数就不是分辨率优度的一个特别恰当的度量，这是因为没有考虑分辨率矩阵的非对角元素是靠近还是远离主对角线，都进行等量的加权。假如存在一种自然顺序关系，那么喜欢比较大的元素都在主对角线附近（如图）。这样分辨矩阵的每一行就是局部化的平均函数。当利用狄利赫莱展布函数计算广义逆时，经常会出现旁瓣，也就是分辨矩阵中大幅值区域远离主对角线。宁愿找到一个没有旁瓣的广义逆，甚至不惜加宽靠近主对角线的非零元素的分布范围，这是因为具有这样分辨率矩阵的解有可能解释成一个在空间上相邻的模型参数的局部化平均值。第二节 地球物理反问题的广义线性逆 (a)(b)旁瓣与 Ba
23、ckus-Gilbert 展布函数 因此，对展布的度量添加上一个加权因子,w i j，根据R的第,i j 个元素离开主对角线的距离对其进行加权。这种加权优先选择呈“尖峰”形，或“类似”形分辨矩阵进行加权。如果自然顺序关系是简单的线性关系，那么选择2,w i jij可能是合理的。如果这样，顺序关系是多维的，就需要选择一个更复杂的加权因子。选择一个使对角元素的权等于零，即,0w i i，并且使,w i j 非负以及使,w i j 关于 i 和 j 对称的展布函数通常是比较方便的。这一新的展布函数通常称为 Backus-Gilbert 展布函数，其具体形式为 第二节 地球物理反问题的广义线性逆 亚定
24、问题的 Backus-Gilbert 广义逆 该问题类似于通过使模型分辨率的狄利赫莱展布取极小来导出最小长度解。因为当该问题是亚定时（因此数据分辨率具有比较小的展布），它非常容易满足数据，所以找出一个仅仅使模型分辨率的展布极小的广义逆。现在要寻找使模型分辨率的 Backus-Gilbert 展布极小的广义逆。因为没有对模型分辨率矩阵的对角元素进行加权，所以还要求所得到的模型分辨矩阵满足方程 这一约束保证分辨矩阵的对角元素是有限的，并且每一行都是对真实模型参数有影响的单位平均函数。第二节 地球物理反问题的广义线性逆 亚定问题的 Backus-Gilbert 广义逆 把分辨矩阵的每一行的展布记为k
25、J，然后代入分辨矩阵的表达式就得到 式中的ijkS定义为 约束方程1ijiiR 的左端也可以用广义逆来表示 式中的ju定义为 第二节 地球物理反问题的广义线性逆 亚定问题的 Backus-Gilbert 广义逆 利用拉格朗日乘子，就能求解使kJ 对广义逆的各元素都取极小问题（在给定的约束条件下）。拉格朗日函数定义为 其中2是拉格朗日乘子。对广义逆的各元素求导数并令其等于零，即 上式必须同原来的约束方程联立求解。把gG的第 k 行作为一个向量，并且把ijkS用下标为 i，j 的矩阵表示，就可以把这些方程写 出矩阵方程，即 这是一个 11MM解方程组，对该方程组进行求解就得到广义逆的第 k 行的
26、M 个元素和一个拉格朗日乘子。第二节 地球物理反问题的广义线性逆 亚定问题的 Backus-Gilbert 广义逆 可以用线性代数中的“镶边法”的一个变种以显式形式来求解该矩阵方程。用这种方法可以构造上述矩阵的逆矩阵。假定该逆矩阵存在，并且可以分为MM阶方阵 A 和一个向量b，以及一个标量c。于是利用子矩阵相等，便可求得未知子矩阵 A，b 和 c。由TijkS AbuI，得1TijkSAIbu,由Au0，得 11TijijkkSSubuu和11TijijkkSSbuuu,由TTijkSc bu0，得 11TijijkkSSubuu和11TijkcSuu。若 A、b 和 c 已知，最终结果为 第二节 地球物理反问题的广义线性逆 协方差大小的计入 可以修改确定 Backus-Gilbert 广义逆优度的度量，使之包含模型参数协方差的大小。利用考虑狄利赫莱展布时所用的那种度量，优度的度量为 01是加权因子，它决定了在度量广义逆的优度时模型分辨率和协方差的相对贡献的大小。第 k 行的优度kJ 为 其中，ijkS由方程1covijijuijkkSSd确定。第二节 地球物理反问题的广义线性逆
