反演实验三

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