线性方程组迭代解法在重力学中的应用

1、线性方程组迭代解法在重力学中的应用肖* 200*9 *学院地球物理学摘要：在重力学中，许多利用地表观测值求内部质量分布和断层产状的反演问题通常存在积分解算。这些积分一般都非常复杂，很难用解析法直接运算，但是我们可以通过线性化来求其近似解。本文讨论的就是重力学中的积分问题通过线性方法解算的内容。关键词：线性代数；迭代解法；重力学一、 重力学中的线性问题1、 重力反演问题在地球重力学中，反演问题是指，由实测重力异常及其倒数的数值大小、空间分布和变化规律，定性和定量的推断地下客观存在的异常地质结构、构造和地质体的形状、产状及剩余密度分布。反演问题的流程可以描绘为：观测数据反演数学物理模型场源模型参数

2、一般地，场源函数与地球表面物理场之间存在着积分关系： （）其中，是场源的作用区域，称为积分核函数（又称格林函数），它刻画“场”与“源”之间的“关系”。若函数为已知函数，积分核函数根据场与场源的物理实际给出，求解场源的物理问题就是反演问题。为了叙述方便，将写成算子形式。若存在逆算子使得，则称m为的（反演问题意义下的、广义的）解。一般情况下，积分算子与微分算子互为逆算子。2、 重力线性问题的线性化求解积分方程的解析求解或精密求解时比较困难的事情，对此，人们往往采用迭代方法或线性近似方法。先考虑线性问题，设为一组归一化的正交多项式基，则一维问题可以变为：（1.2.1）系数待定，带入积分得：(1.2.

3、2)系数为已知。为线性的，称用该方法求解（也就是求解系数）的问题为线性问题。在上述方程中，若已知足够多的观测值，则可能球的所有的。在实际中，不可能使用无限多的观测值去求定无限多个系数，实际中常常使用有线展开方法做的近似，即设关于基函数的选择问题可以这样处理，由于M个观测数值处的坐标已知，则为一组函数，由该组函数构成基函数，记：为常系数，为正交基函数，若则为标准正交基函数。如果在方程中G为一个维向量，则方程演变为一个简单的矩阵方程，即:()方程（）也可以看作是方程（）离散化的结果，这样，反演问题就变成了求解阶线性方程组的问题。3、 重力反演存在的问题反演问题的主要内容有三方面，其一是解得适定问题

4、，包括借的存在性、唯一性及稳定性。其二是反演问题的求解方法。其三是反演问题解的评价。在重力反演中，解的存在性已经被大量事实所证实。然而反演问题的解具有不唯一性，为说明解得非唯一性，让我们来考虑零向量和零空间的概念。若线性反演问题有两个解和，其解非唯一，则有：，两个方程相减，得因为假定这两个解是截然不同的，所以它们的差是非零的，而称矢量为零矢量，而由零矢量组成的空间称为零空间。由此看出，任何具有零矢量的线性反演问题的解都是非唯一的。如果（这里表示特解）是的一个非零解，则也是一个非零解。其中是一个不为零的任意常数。若给定线性问题有q个独立的零矢量，则其一般解为: (1.3.1)总之，反演一组观测数

5、据就是一个，而中的任意一个都可以加到上而仍然拟合观测数据，因而使解变的非唯一。现在回到方程，且设，r为G的秩，显然有以下四种情况：（1） 当M=N=r，且存在时，方程为适定问题，有唯一解。（2） 当M>N=r时，方程为超定方程，无常规意义下的解，但有最小方差解。（3） 当N>M=r时，方程是欠定方程，无常规意义下的解，但有意义下的解。（4） 当min（M,N）>r时，方程的问题是混定问题。这时，只有同时在两种限制条件之下，方程才有解。下面将讨论每一种情况的具体解法。（1） 适定方程对于适定方程，直接使用下文所介绍的线性方程迭代解法求解即可。（2） 超定方程对于超定方程，最简单

6、、最常用的反演方法是最小方差法。(1.3.2)r为G的秩。设e是观测数据g与理论计算值Gm之误差向量，则方差（即目标函数）为 (1.3.3)最小方差解必须满足：所以：(1.3.4)注：当有零特征值存在时，方程是奇异的，无法求解；当的特征值很小时，方程是病态的，会使解变得极不稳定，使反演非常困难，甚至无法收敛。这也是在重力资料反演中经常遇到的问题。（3） 欠定方程从线性代数理论可以知道，此时有无限多个解能满足方程，且其误差均为零。这是因为，虽然数据提供了一些确定参数模型的信息，但其数量不足以全部确定模型参数，或为题够确定模型参数足够充分的信息。因此，解是不唯一的，甚至有无限多能拟合观测数据的解。

7、为了求得重力反演问题的一个解，我们必须从无限多个能拟合观测数据的解中，挑选一个我们所需要的特定解。因此，我们必须加上一些它未包含的信息，这种附加给反演问题的信息叫“先验信息”。在补充先验信息的时候，我们需要讲究“择缺补充”的原则。第一类是补充待求参数的物理性质和可能的数值范围，如地球密度的非负性，且根据地球物理常识可以先定在0.016.0。第二类是先验信息来源于其他一致的地质、重力资料。比如反演地球基底的埋深，油层的厚度，金属矿的属性等。第三类，某些参数比其他参数对解决重力反演问题更重要，此时可以对模型参数经行加权，在一定权系数约束下求解。第四类，假定地球物理模型最简单。这里的最简单是指在保留

8、实际地球物理模型基本特征不变的情况下，最地球物理模型的一种简化。借的长度，比如说借的欧几里得长度为最小的模型，应该是一种简单的模型。当然，这里定义的简单不一定处处非常合理。因此又出来了其他形式的最简单模型，如该模型的变化为最小的范数意义下的模型等。设为一切顶问题，此时的目标函数，在的约束之下有极小，即根据极致理论，必须引入拉格朗日算子将条件极值问题化为无条件极值问题。因此目标函数应为：(1.3.5)显然求上述目标函数的极小值可以化为求：故则 (1.3.6)注：在欠定问题求解中也存在这“奇异”“病态”的问题。此处奇异是指的特征值中有为零的问题；病态问题只是中有小特征值的问题。这都是在解欠定问题时

9、必须认真对待的。（4） 混定方程鉴于混定问题的复杂性，不难想象其目标函数中应该兼有方差项和模型长度项，即(1.3.7)求解可以得到： (1.3.8)此处，称为阻尼系数或加权因子，它决定预测误差项和模型范数长度项再寄笑话目标函数E时的相对重要性，如果阻尼系数足够大，则问题以欠定问题为主，反之则以超定问题为主。二、 线性方程组的迭代解法概述科学研究中及生产实践中，很多问题都可以归结为线性方程组的求解，高效求解线性方程组成了许多科学计算及实践的核心之一。解线性方程组的传统方法是直接法，理论上可以得到方程组的真解，但是很容易出现大数“吃掉”小数、不稳定的问题。线性方程组的迭代解法是一种间接解法，它利用

10、某种迭代过程去逐步逼近准确解，而求出方程组具有指定精确度的近似解的方法。它可以有效地提高求解的速度。迭代法的另一个突出特点是可以充分利用和保持系数矩阵的稀疏性，只需知道系数矩阵与向量乘积的计算法则，而不必知道具体的系数矩阵，这样可以节省内存的开销。可见，用迭代法求解线性方程组，无论在时间上和空间上都有很大的优势。以下，介绍五种常用的迭代方法。1、 雅可比(Jacobi)迭代法对于n阶方程组中的每一个方程（）可以改写为未知向量x的分量的形式：(2.1.1)设方程的准确解为向量，任取一个向量作为的初始近似，将带入（）式的右端，求出的结果记为，称为1次近似。一般的，在求出了第m次近似后，在带入（）式

11、的右端便可得到m+1次近似(2.1.2)给定初值，雅可比迭代公式格式为：(2.1.3)定义1：如果向量序列中的向量每一个分量当时都趋于向量的对应分量，则称是该向量的极限，记为在一定条件下，对任意初始向量，按迭代公式（）求出的向量序列的极限存在且等于方程组的解，这种用迭代格式（）求线性方程组近似解的方法成为雅可比迭代法。2、 高斯赛德尔（Gauss-Seidel）迭代法在雅可比迭代公式中，每算出的一个分量便将其代替该次迭代中以下分量的计算式中的相应分量，便得到另一种迭代公式：（）在给定初始向量，按迭代公式（）求出的向量序列的极限存在且等于方程组的解，这种用迭代格式（）求线性方程组近似解的方法称为

12、高斯赛德尔迭代法。3、 超松弛迭代法（SOR方法）若第次分量迭代已经完成，记第次迭代的第个分量为；第次向量迭代的前个分量已经算出，由高斯赛德尔迭代法算出的第个分量记为：(2.3.1)如果在高斯赛德尔迭代法中引入加权因子，将和进行加权组合作为第m+1次迭代的第i个分量，就得到了松弛迭代法的迭代公式：（）用（）式求方程组解得方法称为带有松弛因子的松弛迭代法：当时称为超松弛迭代方法；时成为低松弛迭代方法；时即为高斯赛德尔迭代方法。4、 对称超松弛迭代法（SSOR方法）如果我们先按自然次序（）用向前的SOR方法逐点计算：(2.4.1)然后再按相反的次序（）用向后的SOR方法逐点计算：(2.4.2)5、

13、 共轭梯度法（CG方法）对于线性方程组的求解问题，定义二次泛函(2.5.1)则可证明求线性方程组的解等价于求二次泛函（）的极小点。由此，给定了初值，按照某一方向去求式（）的极小值点，就得到先一个迭代值，再由出发求最终逼近精确解。若取求最小值的方向为在（）处的负梯度方向，就是所谓的最速下降法。然而理论和实际计算表明这个方法的收敛速度较慢，共轭梯度法则是在处的梯度方向和这一步的修正方向所构成的二维平面内，寻找使减少最快的方向作为下一步的修正方向，即求极小值的方向。计算公式为：(2.5.2)再逐次计算：（表示x,y的内积）可以证明当时，有=0，=0。从而，形成一组共轭向量组；形成以正交向量组。后者说明若没有舍入误差的话，至多n次迭代就可以得到线性方程组的精确解，然而在实际计算中一般都有舍入误差，所以并不是真正相交，所以也只能得到近似解。6.线性迭代的收敛性讨论设使方程组的解，对于给定的初始向量，若由于某种迭代方法产生的向量序列有，则称该方法收敛，否则方法发散。对任意初始向量及任意右端向量f，由迭代产生的迭代向量序列收敛的充要条件是谱半径r(B)<1判别条件1：若|B|<1, 则迭代对任何初始向量都收敛.判别条件2：如果A为严格对角占优阵，则其 Jacobi迭代和Seidel迭代对任何初始向量都收敛。 判别条件3： 如果A为对称正定阵