《地球物理反演概论》上机实验报告实验三：迭代法求解地震层析成像问题
姓名：
学号：
专业：地球物理学
指导教师：邵广周
完成时间：2017.12.21
一、实验内容
利用ART 及SIRT 迭代算法实现下图所示的地震层析成像问题。
⎥⎥⎥⎥
⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=
⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣
⎡020*******
0000
000200020002100100100010010010001001001111000000000111000000000111987
6543
21m m m m m m m m m
二、实验要求
编制相应的程序，在计算机上实现ART 及SIRT 迭代算法。
将ART 及SIRT 算法的反演结果与实验二的结果进行对比分析，比较各反演方法的优缺点。
三、算法原理
对于线性反问题Gm=d ，系统的每一行对应着一个n 维超平面，共m 个超平面。
Kaczmarz 算法的基本原理是：从事先给定的初始模型(0)m 开始，通过将初始模型投影到由G 的第一行定义的超平面上得到(1)m ，然后再将(1)m 投影到由G 的第二行定义的超平面上得到(2)m ，以此类推直到将所有m 个超平面投影完毕。
重复上述迭代过程，直到解成功收敛。
Kaczmarz 算法流程： 1、令(0)0m =
2、for i=0，1，…，m ，()(1)
()
11111
i i i T
i i i T
i i G m d m
m G G G ++++++-=- 3、判断解是否收敛，如不收敛，返回步骤
如果Gm=d 有唯一解，Kaczmarz 算法将收敛于该解。
如果系统有多个解，算法将收敛于与初始模型(0)m 最接近的那个解。
特别地，如果(0)0m =
，我们将会
得到最小长度解。
如果精确解不存在，算法得到的解为最佳近似解。
对于收敛速度问题，当由系统定义的超平面接近正交时，算法收敛速度快。
而当超平面接近平行时，算法收敛速度就会非常慢。
ART 迭代算法：
ART 算法是Kaczmarz 算法的一个修正算法。
它对Kaczmarz 修正算法进行了一个粗略近似，即将G 的第i+1行所有非零元素用1替代。
定义：1i 1j i j q m l ++=
∑
第条射线上的所有单元格（）
为第i+1条射线的近似旅行时，则11q d i i ++-为第i+1
条射线的旅行时预测误差。
ART 算法将Kaczmarz 算法中的修正项()11111i T i i i T
i i G m d G G G +++++-用11
1i i i q d lN +++-替代，其中l 为剖分单元格的尺寸，1i N +为第i+1条射线经过的单元格总数。
因此ART 算法的修正公式可写为
()11
(1)
1()
i 1j i 1j i i i j i i j i j
q d m lN
m m ++++-⎧-+⎪=⎨⎪+⎩第条射线通过第个单元格第条射线不通过第个单元格 该公式可进一步修正为：
()11
(1)
11()
i 1j L i 1j i i i j i i i j i j
d q m lN
m m +++++⎧+-+⎪=⎨⎪+⎩第条射线通过第个单元格第条射线不通过第个单元格 其中1L i +为第i+1条射线的真实长度。
ART 算法流程： 1、令(0)0m =
2、for i=0，1，…，m ，计算第i 条射线经过的单元格总数i N
3、for i=0，1，…，m ，计算第i 条射线的实际长度i L
4、for i=0，1，…m -1；j=1，2，…，n ，计算
()11
(1)
11()
i 1j L i 1j i i i j i i i j i j
d q m lN
m m +++++⎧+-+⎪=⎨⎪+⎩第条射线通过第个单元格第条射线不通过第个单元格
5、判断解是否收敛，如不收敛，令(0)()m m m =返回步骤4进行迭代。
否则返回估计解()m m m =
ART 算法的主要优点：
1、是节省内存，我们只需保存射线经过的单元格的信息，而不需记录每个单元格内各射线的长度。
2、与Kaczmarz 算法相比减少了乘法运算的次数。
缺点：计算精度略逊于Kaczmarz 算法。
SIRT 迭代算法：
SIRT 算法是ART 算法的一个变种，其基本思想是将经过第j 个单元格的所有射线的修正量都计算出来，然后取所有射线修正量的平均值作为模型参数的修正量。
具体算法如下： SIRT 算法流程： 1、令(0)0m =
2、for j=0，1，…，n ，计算第j 个单元格经过的射线总条数j K
3、for i=0，1，…，m ，计算第i 条射线经过的单元格总数i N
4、for i=0，1，…，m ，计算第i 条射线的实际长度i L
5、令0m ∆=
6、for i=0，1，…m -1；j=1，2，…，n ，计算
11
11
i 1j L 0i 1j i i i i j j d q lN m m ++++⎧-+⎪
∆=∆+⎨⎪+⎩
第条射线通过第个单元格第条射线不通过第个单元格 7、for j=1，2，…，n ，令j j j j
m m m K ∆=+
判断解是否收敛，如不收敛，令(0)()m m m =返回步骤5进行迭代。
否则返回当前解。
四、数据及运行结果
图1 输入文件
图2 输出结果
五、实验结论和心得
由输出结果可以看出，ART计算结果与模型真值相符，且均方根误差较小，迭代27次，收敛速度较快；SIRT计算结果与模型真值相符程度较ART差，均方根误差较大，迭代110次，收敛速度较慢。
ART计算结果与实验二Kaczmarz计算结果相差不大，而SIRT效果较差。
ART算法的主要优点是节省内存，与Kaczmarz算法相比减少了乘法运算的次数，因此收敛速度有所提高；缺点是计算精度略逊于Kaczmarz算法。
而SIRT 算法是ART算法的一个变种，目的是提高计算精度。
本次试验中，ART的收敛
速度没有明显提高，并且SIRT的计算精度不升反降，可能是由于网格剖分太少，数据量太小，其优势没有体现出来。
本次试验给定初始模型，通过精度及最大迭代次数控制迭代，在计算机上实现了利用ART及SIRT迭代算法求解地震层析成像问题（已知射线路径及走时数据，求取模型参数即慢度），算法流程清晰明了，因此程序简单易懂。
需要注意
和某一条的一点是，在SIRT算法中，统计了某一个单元格经过的射线总条数K
i
射线经过的单元格总数N
，不要将两者混淆。
i
通过本次实验，对ART及SIRT算法的基本原理及流程有了进一步的理解和认识。