地球物理反演教学课件_第1页
地球物理反演教学课件_第2页
地球物理反演教学课件_第3页
地球物理反演教学课件_第4页
地球物理反演教学课件_第5页
已阅读5页,还剩75页未读 继续免费阅读

下载本文档

版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领

文档简介

时间反复无常,鼓着翅膀飞逝地球物理反演地球物理反演时间反复无常,鼓着翅膀飞逝地球物理反演杨长福E-mail:ycffcy6422sina地球物理反演理论InversetheoryinGeophysics第一章绪论

图2反演问题的形成

反演通常从预测物理场响应开始。怎样从观测值确定未知参数?

图3反演模型的建立以沙滩的理想化模型为例。为使问题简化,将表面看作平面,并将地下划分为若干网格,各网格可能是沙子也可能是金子。“×”是计算重力场点的位置。

首先,我们可以编写以密度为输入、预测重力值为输出的程序,一旦有了这样的工具就可了解不同密度值所能取得的观测重力值,可以假定金子是具有相同大小的矩形块体箱子。我们将箱子移动到不同的位置,即变化箱子的深度和水平位置,看看所得的重力值是否和观测值相一致。编写程序的部分工作是定义将要用到的密度模型。可将沙滩简化为完全水平的表面和由不变密度的沙子所包围的密度可变的一系列矩形块所组成。这些矩形包括金块所有可能埋藏的位置。为模拟设置可能埋宝的矩形块的密度等于金子的密度,其余矩形块的密度等于沙子的密度。图3中“×”是计算重力场点的位置。注意该计算值是预测值而不是观测值。接下来应重点解决的是用程序算出宝物埋藏的位置。假定我们将计算重力的程序插入在一更大程序中,使它具有以下两功能:1)通过小矩形内密度金子和沙子的各种组合产生所有可能的模型2)比较计算值和观测值。若存在的话哪些模型和观测结果一致。模型空间和数据空间:在图3模型中,包含45个参数,即各块的填充物(或金或沙)密度。在数学上用各网格密度将模型表示为一45维向量,如(2.2,2.2,19.3,2.2,…),即该向量就表示了一种模型。而且各矩形块只允许取金子或沙子的密度,因此也可考虑只取0和1的45维向量。那么就有245个这样的模型,将这些模型构成的空间称为模型空间。另一方面数据空间由所有可能预测的数据组成,在本例中有5个观测点,故数据空间为5维。1.1太多的模型

本例中有个模型。若以1千/秒个模型的计算速度,则要花1100年的时间,要完成任何一种可能答案几乎都不可能。1.2无唯一答案本例中有45个调节按钮来组合(各单元一个),而只有5个观测值来匹配,很可能有1个以上能最佳拟合的模型,一旦观测中存在噪声,这种可能性就会变成必然性。对一反演问题总有很多可能的答案,这些答案不能由观测数据来区别。1.3不合理的模型根据外部信息,可以认为财宝等同于6个小矩形金块,也可能埋在一完好无损的箱子里,然而由于暴风雨改变了沙滩的结构,或损坏了箱子使金子散开分布,我们不能绝对肯定任何一种情形,也不能肯定是否被人挖走,根据这些信息,可认为某些模型可能比其它模型更正确,就可给出不同金块数的相对概率,例如,如图4。同样由于倾向于认为箱子完好,故偏爱所有金块都在2×3型矩形块内的模型,而更广范围分布金块的可能性更小,定性地可倾向于某些特定概率的模型,甚至在观测前就如此。如图5所示,这种差别在半定量法很难捕捉。先验(priori)信息和后验(posteriori)信息:独立于观测数据的信息叫先验信息;从观测数据和先验信息所推断的信息叫后验信息。我们所谓的合理信息可等价为独立于观测数据的地下信息,应作为观测资料的附加信息而使用。图4先验条件:预先知道埋藏在沙滩下面的金块情况。如存在被挖走即金块数为0的可能性,而最有可能的是金块数为6左右,但决不会是3等。这些信息表示了独立于重力观测数据的信息,是事先知道的,故称先验信息。

1.4带噪声观测数据大多数观测都受噪声干扰,且重力观测更精密。若有两个模型所产生的预测值是在观测值的误差范围内,就不再强调一个模型可能比另一模型拟合得更好。很清楚要完全了解观测所需要告诉我们的东西,就必须考虑观测噪声。

1.5沙滩不是一个理想模型一个棘手的问题是真实的沙滩不是所要考虑的模型。因为沙滩1)是三维的2)表面不规则3)除沙子和金子外,还有其它物体4)附近有海洋,另外沙滩所寓附的地球除自己的质量外还受月球和太阳的吸引。5)其它等等这些效应使得需要拟合的观测资料需要校正,这有可能增加观测误差。1.6小结反演理论是根据观测数据(通常是由遥感得到)得到的物理方程作出推断的问题。既然所有数据都有不确定性,这些推断就具有统计性。而且所记录的数据仅仅是有限的带噪数据,且物理机制是由连续方程来模拟,故只要有一个模型能拟合观测数据,那么就会有无限多这样的模型。为使这些问题定量化,须回答三个问题:1)观测数据精度如何?即“拟合”数据是什么意思。2)正演的精度如何?即是否考虑了所有对观测数据有较大贡献的物理因素。3)对独立于观测数据的机理知道些什么?由于对一物理模型的任何足够好的参数化,都可能存在也能拟合数据的不合理模型,就必须存在舍弃这些不合理模型的系统方法。1.7例子下面列举一沙滩计算的例子,如图6所示,随意调整金子/沙子矩形格子,以图计算的响应拟合观测数据。对这种特殊的计算,真实模型有6块金砖(如图6),图6b)显示一也拟合观测数据的模型,尽管预测数据和观测数据间的差别不是严格为0,但在观测中有噪声,故可视为拟合得足够好。因此看到两不同模型几乎拟合观测数据一样好。图6

a)金块的真实分布

b)能拟合观测资料的不合理模型

1.7不是反演的简单反演问题现在看看另一反演问题:从物体的重量和体积估算其密度。尽管这对反演者似乎是一个太简单而无任何兴趣的问题,但将向你展示它是和根据地震观测模拟地球三维弹性结构完全一样的令人头疼的理论问题。图7是一不规则物体K,我们的目的是估算其密度,密度为标量(scalar)如7.34,用ρ表示图7物理性质未知的物体图8未知块体的体积测量图9未知块体的质量测量由于可能或不能直接用天平测出质量,在此情形下上,常测量作用在该质量上的重力,然后求出质量。为测出密度,先要测出体积和质量,如上图所示,然后计算出密度。设测出结果是:则可得:在大多数情况下,测量就到此为止,但并不知道测量误差。1.8质量测量误差(噪声)为简单起见,设体积测量无误差,而集中考虑质量误差。经大量测量,测量结果的密度分布如图所示。真实密度为5.2时的测量密度概率分布

1.9如何得到测量结果?设观测密度为,为真实值,表示真实密度为时测得密度值为的条件概率。首先要记住事先并不知道真实的密度值,但是若反复测量,就可给出概率分布图。其次还要注意上图中真实的密度值并不严格落在观测概率分布的峰值上。这可能是测量中某种系统误差造成的,譬如标度偏差、A/D转换器造成的误差、甚至定向的微风等。需要区别理论误差(模拟误差)和随机误差。例如测量铜的密度,理论上总是假设样品完全由铜组成的,但实际上它是不可能的,总会含有一些其它杂质,由此而引起的误差是理论误差的一个例子。在此情形下理论上会有些错误。而随机噪声(又称随机误差)是在反复测量时所发生的微小扰动。不幸的是在实际中,很难区分这两种误差的差别,几乎没有可重复的实验。所以在重复一实验时,实际上总是将一些小的变化带到所作的假设中去了。如重新拿起一样品放在天平上称其质量时,可能有不同的灰尘沾在上面,或将其放上天平时,可能触到天平而轻微地破坏了天平的平衡。就目前而言,理论误差和随机误差的判断有点随意和随机而定,至于哪些误差是我们感兴趣的,哪些不是这里不详细讨论,请阅读文献Whatis

Noise?(J.A.ScalesandR.Snieder.Whatisnoise?Geophysics,63:1122-1124,1998.1.10正演和反演(forwardandinversion)减是加的反问题,除是乘的反问题,开方是乘方的反问题。如给出答案6,那么提出的问题是3+3还是15÷5或?由此可见反问题有无穷多解,这是反问题的内在特性。正演:由已知模型,计算模型的观测响应。数学描述:已知f(m),求y,即计算y=f(m)。m—模型参数,f—模型和响应的函数关系,y—计算的模型响应(数据)反演:由已知观测数据(响应),求模型。数学描述:已知y,求m,即m=f-1(y)关系:正演是反演的前提和基础,反演是正演的逆过程。1.11如何解反问题?Ingeneralworkwelookforanewlawbythefollowingprocess:Firstweguess

it.Thenwecomputetheconsequencesoftheguesstoseewhatwouldbeimplied

ifthislawthatweguessedisright.Thenwecomparetheresultofthe

computationtonature,withexperimentorexperience,compareitdirectlywith

observation,toseeifitisworks.Ifitdisagreeswithexperimentitiswrong.In

thatsimplestatementisthekeytoscience.

—RichardP.Feynman猜测—计算—比较修改猜、算、比、改2.反演理论研究内容a)解的存在性:b)模型构置:c)解的稳定性和非唯一性:d)解的评价:3反演史和应用X-raytomography(XRT)(CT=computedtomography:切片图).PositronEmissionTomography(PET)(正电子放射层析成像)Wholeearthtomography:Localearthquake®ionaltomography:用地震波对地球患者的CT成像4、反演问题分类线性非线性线性化目前非线性问题的反演主要是建立在线性化基础上,即非线性问题的线性化反演,但是一种近似。线性化是一种近似(无奈之举)参数化是一种选择(建立离散化模型)非唯一性是一种本能原因:问题本身性质误差观测信息有限现有正演理论及计算、反演方法不完善第一章地球物理反演概述一.常用的基本数学概念1.矢量(vector)又称向量a)定义:在物理上曾定义为有大小和方向的物理量(狭义)。现给出广义上的矢量概念:由若干个参数所描述的数学模型或物理量,常记为:,。

—第i个向量分量,N—向量维数,—N维向量。在算法中称N维数组。在数学里常视为矩阵的行或列,又称行矢量或列矢量,故除它自身的特点外,还具有矩阵的特点。事实上矢量在数学中,完全可视为只有一行或一列的矩阵,同时矩阵也可视为若干相同大小(维)的矢量组合,因此矢量满足矩阵的运算规则和书写规则。需要注意的是,这也是书写造成混乱的原因。和矢量相对的是标量(scalar),即无方向的物理量,或仅用数量表示的量,在数学上表示为单个数据。

方向:用的单位矢量(unitvector)表示:b)矢量大小(长度、模)c)线性矢量空间(LinearVectorSpaces)所有满足线性运算法则的N维矢量所构成的空间d)矢量运算两个维矢量的加法:若上式等于0,则称向量和正交,即垂直。

标量乘法(内积,innerproduct,投影,Projecting):e)矢量的线性相关和独立(LinearDependenceandIndependence)对n个同样维数的矢量:若存在不全为0的系数使成立,则此n个矢量线性相关,反之只有全为0的情况下,上式才成立,则称n个矢量线性无关(独立)。即能被别的矢量表示的矢量叫线性相关,否则线性无关。线性方程有唯一解的充要条件就是矩阵组成的所有行(列)矢量线性无关。n维线性矢量空间最多只能有n个线性无关的矢量,任何个数超过矢量维数n的矢量必然线性相关。f)关于共轭(正交)的向量:设向量、、…、是个非零向量,为对称正定矩阵,若对任意有:则称向量、、…、是关于共轭(正交)的向量组。共轭向量线性无关。2矩阵(Matrix)a)定义:将一维空间向量变换为维空间向量的线性变换(或映射),叫线性变换的系数矩阵,简称矩阵,记为:当时称方阵。

若,则矩阵乘法:注意:对角矩阵:单位矩阵:矩阵转置(transpose):行换成列,或列换成行。如矩阵,其转置矩阵记为:两矩阵乘积的转置:b)特征值和特征向量(EigenvaluesandEigenvectors)矢量经矩阵变换后,仅是增大或缩小,方向不变:则称的特征值,称的特征矢量。由得特征方程:

称特征多项式

例特征多项式为:求得特征值:代入特征方程可得:和可求得两特征矢量:注意:特征矢量不是唯一的,因为可任意乘一因子而不改变其特征矢量的性质。c)矩阵的秩(rank):矩阵中不为0的代数余子式的最高阶数称为矩阵的秩。也即是矩阵中所包含的线性独立特征矢量的最大个数。当矩阵的秩为,则有个正特征值,其余特征值为0d)正交矩阵和正交变换:若线性变换保持模不变,即有:,则线性变换叫正交变换,相应的矩阵叫正交矩阵。为正交矩阵的充要条件是:

AT=A-1

e)正定矩阵:对任一组不全为0的x1,x2,…,xn,二次型:f(x1,x2,…,xn)的值总是正的,则其所对应的对称矩阵,即

:f(x1,x2,…,xn)=xTAx中的叫正定矩阵,其特征值均大于0。f)逆矩阵方阵A的逆矩阵A-1(普通逆):AA-1=A-1A=IA-1=其中A*为A的代数余子式组成的矩阵,若︱A︱不等于0,则称A为满秩矩阵或非奇异矩阵。对任意满秩矩阵(行满秩、列满秩和方阵)存在左逆和右逆:左逆(leftinverse):

右逆(rightinverse):

行满秩矩阵,存在右逆,列满秩矩阵,存在左逆。只有当矩阵为方阵时,才可能同时存在左逆和右逆,且等于普通逆。须指出:满秩矩阵的逆(无论是左逆还是右逆)都不是唯一的。对满秩矩阵,其左逆和右逆的一般形式是:左逆=

右逆=

式中V为适当阶数的任意矩阵。例若选矩阵V为单位矩阵,则:左逆=

=若选择,则可得另一左逆==经验证,这两左逆都满足左逆定义公式,说明了逆矩阵的不唯一性。满秩分解:对m×n阶非满秩矩阵A,总存在它的一个满秩分解,即存在一m×r阶行满矩阵C和一r×n阶列满矩阵D,使广义逆(generalizedinverse):设对满足以下四个条件:则称为的Moore-Penrose广义逆,记为以示和普通逆的区别。广义逆具有唯一性。

定理(Theorem)1):对一n个方程m(n≤m)个未知数组成的方程组:Ax=y,A为n×m阶矩阵,当且仅当A的秩r=n时,存在一m×n阶右逆矩阵C,使定理(Theorem)2):对一n个方程m(n≥m)个未知数组成的方程组:Ax=y,A为n×m阶矩阵,当且仅当A的秩r=m时,存在一m×n阶右逆矩阵B,使定理3)(满秩矩阵分解):设A为一n×n阶矩阵,S为A的n个线性独立特征矢量组成的列矢量矩阵,则A可对称化为其特征值组成的对角矩阵,且=定理4)对一n×n阶矩阵A,若存在相应于n个不同特征值的n个特征矢量,则这些特征矢量线性独立。定理5)

实对称矩阵A可分解为由正交矢量Q组成的矩阵和实特征值组成的对角矩阵的乘积(正交分解):3范数(norm)对向量或由向量所组成的矩阵按某种规范所定义的数,分向量范数(又称L范数,可理解为向量大小、长度等)和矩阵范数。向量范数(范数):常用的有三种范数,分别称、和范数。向量椭球范数:在空间中,定义为:矩阵范数(F范数)1)2)3)方阵的条件数(又称状态数,conditionalnumber)定义方阵A的条件数:注意它和范数的定义有关,通常用。记为对方程的性能可由A的条件数来表示:、分别为的最大和最小特征值,状态数越小,方程性能越好,反之方程性能越差,当接近无穷大时,称为病态(ill-posed)方程。方阵的特征值有为0的问题称奇异问题,当特征值中出现小特征值时为病态问题。4方程解的不稳定性:

1.11x1+x2=3.11(1)x1+0.90x2=2.80其解为x1=1,x2=2,但当右边存在误差(-0.01,0.01)时有1.11x1+x2=3.10(2)

x1+0.90x2=2.81则解为x1=20,x2=-19.1。这种因观测较小误差而引起解较大误差的性质称解的不稳定性。b)吉洪诺夫(TihonovA.H)正则化(regularization)单一密度界面的反演公式:D是一经典的解不稳定问题。正则化的目的就是使不稳定问题得到满足某种条件(先验的或约束的)下的稳定解。换句话说正则化就是对方程加上某种条件(正则化条件),使其得到该条件的稳定解,该过程也叫正则化过程。例:对Ax=b其最小二乘解为:最小长度(范数)解为:阻尼最小二乘解:其中为阻尼因子,为I单位矩阵。下表为上例中方程(2)在取不同值的阻尼最小二乘解k12345610.10.010.0010.00010.00001x11.23981.50921.54311.55071.59201.9920x21.11641.35891.38871.38721.34200.89790.82990.10190.01720.01380.01370.0134e1.66830.36250.04510.00770.06120.059710.91560.81870.81780.82390.88511.4828表中为观测值与理论值的均方差,e为解(最小长度解)与先验估计的均方差,为反演解与真实解(1,2)T的均方差。

5多元函数的梯度、海森矩阵、雅可比矩阵设函数对变量的一、二阶导数存在,则梯度:海森(Hessen)矩阵(对称):雅可比(Jacobi)矩阵:设向量函数为:则称雅可比(Jacobi)矩阵。

6多元函数的泰勒展开借助于一元函数的一、二阶导数,可导出元函数的泰勒展开。设n元函数的泰勒展开:略去高阶小量,得一次近似:和二次近似:7解Ax=b(A为对称正定)与求二次函数与

的极小值等价。

求或是Ax=b的解的充要条件是:是的极小点。证明:必要性:由函数极小的必要条件必有:即是Ax=b的解。反之(充分性):即是的极小点。8线性(linearity)和非线性(non-linearity)当函数满足以下两个条件,称为线性函数或线性泛函。反之称非线性函数。二、线性反演的一般理论1正演和反演(forwardandinversion)的一般概念正演:由已知模型,计算模型的观测响应。数学描述:已知f(m),求y,即计算y=f(m)。

m—模型参数,f—模型和响应的函数关系,y—计算的模型响应(数据)

反演:由已知观测数据(响应),求模型。数学描述:已知y,求m,即m=f-1(y)关系:正演是反演的前提和基础,反演是正演的逆过程。2正反演中另一更直观的定义:正:(1)

为数据空间中的一点(即一组观测数据),为模型空间中的一点(即用一向量表示的一个模型,向量分量为模型参数),为泛函算子,反映由模型空间到数据空间的映射,可以是积分算子或微分算子,也可是矩阵或函数,当它是线性算子,则是线性反演问题,反之若它是非线性算子,则是非线性反演问题。该方程的解就是反问题:模型空间模型空间和数据空间的映射关系图数据空间3.地球物理反演问题地球物理反演问题是地球物理学中的一个核心问题。场和场源及两者的关系是地球物理学中的三个基本要素。场和场源的关系可用三种基本形式表示,除上述矩阵方程形式(1)的表达外,还有微分方程形式、积分方程形式。(1)还可写成:积分方程形式:(2)为场的空间分布函数(即观测数据),为场源的模型函数,为积分核函数,又称Green函数,表征场和场源的关系,反演问题就是确定待定的模型函数,反演问题就是求解积分方程问题。为与场源性质及分布有关函数,为微分算子,对不同的问题代表不同的算子。为场的分布函数。微分方程的正演问题是通过给定的边界条件或初始条件,确定方程的定解,以求得那些不可直接测得的场分布或重建方程系数,进而确定模型参数。地球物理反演问题一般都是非线性问题。微分方程形式:4非线性问题的线性化

方法:利用多元函数的泰勒展开可将非线性问题转化为线性问题。设一物理(地质)模型由个模型参数确定,则第j个观测数据是该n个模型参数的函数,即可表示为:由泰勒定理(略去非线性项):即(3)若观测数据为m个,则观测数据可表示为m维向量:其中分别表示不同的观测数据。

同时有m个类似(3)的方程,这m个方程写成向量,就是:即:

(4)

其中称雅可比(Jacobi)矩阵,又称灵敏度矩阵(Sensitivity)。因此,非线性反演问题,经泰勒展开,可转化为线性矩阵方程组。另一线性化途径是从目标函数(方差函数)出发,经泰勒展开,也可线性化.。目标函数:设分别表示不同的观测数据,则观测数据可表示为向量:并令

为各观测数据所对应的“理论观测场值”(通常称理论计算场值),并写成向量:则所有观测值和理论值之间的差的平方和称为目标函数:

(5)

将它在展开并略去高次项即(非线性项),得:写成向量形式:(6)即海森矩阵。

,其中对(6)取极值,由极值条件:

可得:

(7)形式完全相似于(6)。连续模型的离散化和线性化:式中写成矩阵方程有:(8)其中,反演之父R.Parker(1970),将反演归为四个基本问题:

1)解的存在性:解是否存在2)模型构置:即建立模型求得稳定解。3)解的唯一性:解是否唯一。

4)解的评价:如何从非唯一解中选出真实信息(或正确解)。第二章线性方程求解方法

一、方程分类设有M个方程(即有M个观测数据),N个未知数(即N个模型参数),则系数矩阵为M×N阶,其秩为r.1)若M=N=r,适定(good-posed)方程,信息不多不少,方程有唯一解。2)若M>N=r,超定(overdetermined)方程,观测资料存在多余信息,方程无常规解,但有最小方差解。3)若N>M=r,欠定(overdetermined)方程,观测资料提供信息不足,方程无常规解,但有最小长度解。4)若min(M,N)>r,混定方程(奇异方程:出现特征值为0;和病态方程(ill-condiontionedorill-posed):出现特征值很小),观测资料足够多,但仍不能提供足够的独立信息,方程无常规解,只有最小方差和最小长度共同约束下的解。方程有常规解,则称相容方程,若无常规解则称不相容的或矛盾方程。二、适定和超定方程的求解方法

设有适定(welldetermined)和超定(overdetermined)线性方程d=Am,按常规方法求出精确的解极不稳定,毫无意义。通常可求观测数据和预测数据误差平方和最小(一种正则化手段)时所对应的解(即模型参数)。为此使方差最小:的解称最小方差解或最小二乘解(又称最小L2解)。由(2-1)展开得:(2-1)由极值条件:

(2-2)即相当于在(2-1)两端乘,即使(2-1)正则化。得最小方差解或最小二乘(LeastSquares)解:注意:当(N×N阶)的秩为r<N,即为奇异(非满秩,有0特征值)方阵,则方程是奇异的,此时(2-2)无常规解,即无法求解,也称发散。而当的特征值很小(接近0)时,即为病态方程,(2-2)的解极不稳定。即是说,有些超定方程是不能用这种方法求解的。三、欠定(underdetermined)方程的求解方法在欠定情况下,即缺少信息。不能用上述最小方差法求唯一解。事实上,对欠定方程,不止一个误差为0的解,即存在无穷多个方差最小的解,因此不能提供足够的信息来唯一地确定方程的解(即模型参数)。为得到这无穷多个误差为0的解中某个有用的解,就必须增加未包括在方程中的信息,称先验信息。增加先验信息原则是:缺什么,补什么,且所补信息一定要正确可靠。先验信息分类:1)模型参数的定义域约束,即未知数取值范围规定。如密度、电阻率、速度必须大于0。若还知更准确的取值范围,还可进一步具体地约束取值范围。2)已知信息(资料),如钻井、地质资料等。3)对参数进行加权,即强调一些参数比另一些参数的重要性。4)“模型最简单原则”:“最简单”在反演理论中用解的“长度”来度量,一般当解的欧几里德长度即L2范数度量下最小,被定义为“最简单”模型。即求在约束下使L2=极小的条件极值问题。因此,目标函数为:由极值条件:可得:代入中可得,则得,所以得最小模型解或最小长度解:(2-4)由此可知是在两边乘以(正则化)得到的方程所求得的解。对比超定方程的解,那里是要求(N×N阶)的逆,这里是求(M×M阶)的逆,同样也存在“奇异”和“病态”两大问题。三、混定(mixeddetermined)方程的求解方法1

阻尼最小二乘法(DampedLeastSquares),又称马奎特(Marquardt)法大多数方程表现为混定方程,既不完全超定,也不完全欠定,或既具超定性质,又具欠定性质,这样不能用最小二乘法也不能用最小长度法求解,但可用两者的结合方法。取两目标函数的线性组合得:(2-5)由极值条件:得:

为加权系数,表示预测误差和模型长度在极小过程中的相对重要性,称阻尼因子或加权因子。若足够大,则第一项可忽略,即使解的欠定部分达到最小,若取0,则是误差达到最小。当取适当值,则得一折衷解:(2-6)显然是最小二乘正则化方程在对角线元素上加一阻尼因子。故这一方法称阻尼最小二乘法,又称马奎特(Marquardt)法、脊(岭)回归法(ridgeregression),其解称阻尼最小二乘解。另一有效方法是广义逆法。2.广义逆法(1)广义逆定义:则称为的广义逆,并记为(2)广义逆性质:广义逆的计算方法:

1)若A满秩,则2)若A为对角矩阵,则,其中3)若A为行满矩阵,则4)若A为列满矩阵,则5)若A为非满秩矩阵,则A可分解为一个m×r阶行满矩阵C和一个r×n阶列满矩阵D,使A=CD,且6)最有效的求广义逆方法:奇异值分解法(3)奇异值:对任意M×N阶矩阵A,设ATA或AAT有r个大于0的特征值,则,,称为A的r个奇异值,且。(4)Penrose奇异值分解定理:设A为M×N阶任意矩阵,且秩为r,则必然存在一个M×M阶正交矩阵U(由AAT的特征向量组成)和N×N阶正交矩阵V(由ATA的特征向量组成),使或其中为M×N阶对角矩阵:,,,为A的r个奇异值,是非增序列。这是奇异值分解的重要特征。因U和V是正交矩阵,所以有:,,故A的广义逆矩阵为:,(2-8)则方程的解为:(2-9)(5)Lanczos奇异值分解定理:设A为M×N阶任意矩阵,且秩为r,则可分解为式中Up为AAT的P(P≤r)个最大特征值对应的特征向量组成的M×P阶半正交矩阵,Vp为ATA的P(P≤r)个最大特征值对应的特征向量组成的N×P阶半正交矩阵,为对角矩阵,其元素为由大到小依次排列的的P个最大特征值的平方根。由Lanczos奇异值分解得到的广义逆为:(2-10)则方程d=Am的的解为:其中Lanczos奇异值分解是Penrose奇异值分解的非0部分:U0、V0是和0特征值对应的特征向量所组成的矩阵,Up、Vp是U、V的一部分,不是正交矩阵,而是半正交矩阵,即:可以证明,由奇异值分解得到的广义逆满足Moore-Penrose广义逆的四个基本条件,且唯一和稳定。当矩阵存在一扰动E,则奇异值变化不超过ETE的最大特征值的平方根。所以奇异值分解是相当稳定的算法。线性方程d=Am的条件数为:(6)奇异值和条件数:最小二乘方程的条件数:很显然,最小二乘正则化过程使方程的条件数增大,加剧了方程的病态程度。阻尼最小二乘方程的条件数为:适当选择阻尼因子,就使条件数相对最小二乘法减小,改善了方程的性能。(7)阻尼最小二乘法(马奎特法)的奇异值分解法(修改后的广义逆法):将阻尼最小二乘方程的解,用奇异值分解结果写出:(2-13)其中为对角矩阵。可见只要将奇异值分解中的对角元素换成就可实现阻尼最小二乘法的奇异值分解。由上式可知,在中的奇异值很小即时,若不加阻尼即则出现分母趋近0的情况,使方程的解极不稳定。(8)小奇异值的截断

由上述分析可知,奇异值越大,对重建观测数据(又称预测数据)的贡献越大,反之越小,当大、小奇异值相差悬殊时,小奇异值对重建的观测数据几乎毫无作用,甚至去掉也不会影响观测数据的重建精度。另一面小奇异值对方程性能影响极大,它越小,对所构置的模型参数越大,即它的一较小的变化会引起模型参数的巨大变化,使解极不稳定。因此要截断小奇异值,这首先由Wiggins(1972)提出,并取得了和阻尼最小二乘法相似的效果。而Lanczos分解定理为截断小奇异值提供了理

温馨提示

  • 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
  • 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
  • 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
  • 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
  • 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
  • 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
  • 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。

评论

0/150

提交评论