反演实验四

《地球物理反演概论》上机实验报告实验四：曲线拟合问题的共轭梯度法
姓名：
学号：
专业：地球物理学
指导教师：邵广周
完成时间：2017.12.26
一、实验内容
利用共轭梯度法实现下图所示的地震层析成像问题。

⎥⎥⎥⎥
⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣
⎡020*******
0000
000200020002100100100010010010001001001111000000000111000000000111987
6543
21m m m m m m m m m
二、实验要求
编制相应的程序，在计算机上实现共轭梯度算法。

三、算法原理
考虑二次最优化问题： 其中，A 为n n ⨯阶的对称正定矩阵，要求A 正定的目的是保证目标函数()X φ收敛且有唯一极小值。

我们可以通过计算目标函数的导数并令其等于零来求极小值，即
()b AX X -=∇φ
极小点处的X 满足：
0=-b AX 或b AX =
因此，求方程b AX =的解等效于求()X φ的极小值问题。

共轭梯度法解最优化问题是通过构造n 维向量基110,,,-n P P P 来实现的，即
0=j T i AP P j i ≠
具有上述性质的向量则称它们是关于矩阵 A 相互共轭的向量。

X 可用向量基展开为如下形式：
()X b AX X X T T
-=
2
1min
φ
∑-==1
n i i i P X α
因此
()⎪⎭
⎫
⎝⎛-⎪⎭⎫ ⎝⎛⎪⎭⎫ ⎝⎛=∑∑∑-=-=-=10101021n i i i T n i i i T
n i i i P b P A P X αααφ
上式可写为：
()⎪⎭
⎫ ⎝⎛-=∑∑∑-=-=-=1
0101021n i i i T n i n j i T i j i P b AP P X αααφ 由于向量关于A 相互正交，上式可简化为：
()⎪⎭
⎫ ⎝⎛-=∑∑-=-=1
010221n i i i T n i i T i i P b AP P X ααφ 上式表明()X ϕ由n 项组成，且每一项彼此独立。

因此只要保证第i 项的系数
i α使该项最小，从而使各项之和达到最小，第i 项为：
i T i i T i i P b AP P αα22-
上式关于i α求导，并令导数等于零，可得使第i 项最小的最优系数i α，即
i
T
i i
T i AP P P b =α 因此，只要我们知道关于A 共轭的一组向量基，则()X φ的最优化问题就非常容易。

那么，如何构造一组共轭向量呢？
共轭梯度算法实际上是通过迭代生成一系列解向量i X ，残差量i
i AX b r -=和共轭向量基i P 。

算法从00=X ，00=r ，00r P =，0
00
00AP P r r T T =α开始迭代。

假设前k 次迭代已得到解向量k X X X ,,,10 ，残差向量k r r r ,,,10 ，向量基k P P P ,,,10 和最优系数k 1
0ααα，，， 。

并假设这1+k 个向量i P 关于A 共轭，向量i r 相互正交，且0=j T i P r j i ≠
令
k
k k k k k k k AP r r P X X αα-=+=++11
上式正好是残差向量的修正公式，这是因为：
()()k k k k k k k k k k k AP r AP AX b P X A b AX b r ααα-=--=+-=-=++11
令k
k k k k
T
k k T k k P r P r r r r 1111
11+++++++==ββ 下面证明k T k k T r r P b =：
下面证明当时与相互正交。

对于任意的，有
由于与之前的所有相互正交，所以
此处由于A 为对称矩阵，所以。

又因，所以有
而和与互为共轭向量，所以 当时我们同样可以得到：
接下来证明当时与相互正交：
()()()
k
T k k T k k k T k k T k k k T k k T k k k k T k k
T k k T k k
T
k k k T r r r r P P A P P r r r AX P P r r AX P P r P AX r P b =++=++++=++=+=+=----00110011ααββ k i ≤1+k r i r k i <()k T
i k i T k i T
k k k i T k AP r r r r AP r r r αα-=-=+1k r i r k T i k i T k AP r r r α-=+01k T i i T k AP r Ar P =1-+=i i i i P r P β()k T
i i i k i T k AP P P r r 110-+--=βαi P 1-i P k P 01=+i T k r r k i =01=+i T k r r ()()0
0111=-=+-=+-=--=-=--+k T k k T k k k k T T k k T k k k k T k k k T k k
T
k k k k k T k k
T
k k k k T k r r r r P b r AP P AP P r r AP P P r r r AP r r r βαβααβααk i ≤1+k r i P
最后证明当时，对于
，有 对于，有
由以上推导可知，算法通过迭代生成一系列共轭向量。

理论上通过次迭代可找到系统的极小解。

而实际上由于舍入误差的存在只能得到一个极小解附近的一个近似解。

共轭梯度算法流程：对于给定的正定对称系统b AX =，给定初始值0X ，令
0,010==-p β，0,00=-=k AX b r 重复如下步骤直到解收敛。

()
()0
0001111
1111111=-=-=-=+=+=+=----+-++-++k T i k i i T i i i T
k k k i i T k i i T k i i T k i i i T k i T k AP P P r P AP r P r P r r r P r r P r αββαββββk i ≤01=+i T k AP P k i <()()()
110
1111111111=-=⎪
⎪⎭⎫ ⎝⎛-=+=+=+=++++++++++i T k i T k i
i i i T k i T k i
T
k k i T k i
T
k k k i T k r r r r r r r AP r AP P AP r AP P r AP P ααββk i =()()()()
()
()
01
111
1
111
11
11111
1111111=-+=
-+=-+=⎪⎪⎭
⎫ ⎝⎛-+=+++++++-++++-+++++k T k k k k T
k k
k T
k k T k k k k T
k k k k T k k T k k k k k
k k k T k k k k T k r r r
r
r r r P r r r r r P r r r P r AP P ββαβββαββααβn
上述方法只适用于正定系统，对于任意系统d Gm =，可通过如下变换，使其
变为正定对称系统。

m G Gm G T T =， 如果令k T k k k s G r Gm d s =-=,，则：
最小二乘共轭梯度算法流程：对于给定的正定对称系统d Gm =，令k=0,
00=m ,00001,0,0,0s G r s p T ====-β，重复如下步骤直到解收敛。

1.如果k>0，令1
1--=k T k k
T k k r r r r β
四、数据及运行结果
图1 输入文件
图2 输出结果
五、实验结论和心得
由输出结果可以看出，共轭梯度法计算结果与模型真值相符，且均方根误差较小，迭代9次，收敛速度较快。

本次试验给定初始模型，通过精度及最大迭代次数控制迭代，在计算机上实现了利用共轭梯度算法求解地震层析成像问题（已知射线路径及走时数据，求取模型参数即慢度），算法流程清晰明了，因此程序简单易懂。

需要注意的一点是，共轭梯度算法要求正定对称系统b AX =，对于任意的n m ⨯方程组n m ≠的一般情况，等式两边需同时左乘A T 即d b A AX A T T ==，转换为n n ⨯的对称正定系统，再进行计算。

通过本次实验，对共轭梯度算法的基本原理及流程有了进一步的理解和认识。