地震走时层析成像方法中大型稀疏矩阵方程解算方法比较.docx
《地震走时层析成像方法中大型稀疏矩阵方程解算方法比较.docx》由会员分享，可在线阅读，更多相关《地震走时层析成像方法中大型稀疏矩阵方程解算方法比较.docx（7页珍藏版）》请在冰点文库上搜索。
地震走时层析成像方法中大型稀疏矩阵方程解算方法比较
地震走时层析成像方法中大型稀疏矩阵方程解算方法比较
X
杨海波
（成都理工大学信息工程学院,四川成都 610059
摘 要:
在地震走时层析成像方法中,对大型稀疏矩阵方程解算方的稳定性和计算效率的研究是一个重要问题。
本文列举了五种用于解决此类问题的常用算法:
SVD、SIRT、Jacobi-ICCG、SSOR-ICCD和LSQR,并在求解一般问题和病态问题的试验中,分析各方法的误差、计算效率和稳定性;进而使用各算法对同一地质模型进行速度成像,并对结果进行比较,给出针对各种方法的客观评价。
关键词:
地震走时层析成像;稳定性;计算效率;误差分析
中图分类号:
P631.4+43 文献标识码:
A 文章编号:
1006—7981（201010—0051—021 问题的提出
地震走时层析成像是根据地震波走时观测数据求解地球内部结构的方法。
在这项研究中,根据对问题的地质了解假定一个初始模型速度模型,用射线追踪方法得到射线的理论走时,解决模型正问题,再借助于模型参数扰动,构成初始模型的理论值和观测值之间的线性方程,反演速度结构,修改初始模型重复上述过程,直到获得满意的结果。
一般而言,地震层析成像反演问题可以归结为解不适定方程组:
Ax=B（1通过构造正规方程:
ATAx=AT
b
（2解出超定方程的最小二乘解,就可以得到这个断面中各个单元的慢度,从而推测出断面中不同物质的分布状况。
目前,解决这个问题主要有三类方法:
¹SVD（奇异值分解法。
ºSIRT（联立迭代图象重建法。
»投影方法。
其中又有CG（共轭梯度法和LSQR（最小二乘共轭梯度法。
在CG法中,广泛运用于实际的是PCG（预测共轭梯度法中的几种算法,本题选用ICCG（不完全Cholisky分解共轭梯度。
本文将用上述方法求解一般问题和病态问题,分析各方法的误差、计算效率和稳定性;进而使用它们对同一地质模型进行速度成像,并对结果进行比较,给出针对各种方法的客观评价。
2 一般问题求解试验
先用规模较小的矩阵方程对各算法程序进行试运行,其目的是检验算法在理论上的正确性。
为此构造一个系数矩阵A的阶数n=5的小型矩阵方程Ax=b（按方程（2,A采用对称正定矩阵的形式:
2.4 2.5 6.9 4.3 1.32.5 5.6 6.1 2.8 7.76.9 6.1 4.9 8.1 5.74.3 2.8 8.1 5.9 1.81.3 7.7 5.7 1.8
3.5x1x2x3x4x5=
4.6
5.1
6.52.43.3该方程的准确解为:
X=（-4.635068,0.279888,-0.670548,4.
287092,0.935946T,各个方法算出的解为x,设d=Ax-b定义了解出解与准确解之间的差向量。
表1列出了准确解和SVD、SIRT、Jacobi-ICCG、SSOR-ICCD、LSQR五种方法得到的解和误差,迭代计算在解的残差小于1.0×10E-6时停止计算。
表1
五种方法解算结果对比
准确解
SVD
SIRT（迭代次数:
55717
Jacobi-ICCG（迭代次数:
5
解
残差解残差解残差-4.635068-4.6350680.000000-4.635068-0.000000-4.635064-0.0000010.2798880.2798880.0000000.279888-0.0000000.2798880.000000-0.670548-0.670548-0.000000-0.670548-0.000000-0.670548-0.0000004.2870924.287092-0.0000004.287092-0.0000004.2870890.0000010.935946
0.935946
0.000000
0.935946
0.000000
0.935946
0.000000
准确解
SSOR-ICCG（迭代次数:
5,
松弛因子:
X=1.0LSQR（迭代次数:
5解
残差解残差-4.635068-4.6350680.000000-4.6350680.0000000.2798880.279888-0.0000000.279888-0.000000-0.670548-0.670548-0.000000-0.670548-0.0000004.2870924.287092-0.0000004.287092-0.0000000.935946
0.935946
0.000000
0.935946
0.000000
从表1可以看出,所有的五个算法能得到非常
接近准确值的解,这说明所有的算法在理论上都是正确的。
但很明显地,五个算法在得到相近结果的同时,所付出的计算上的代价并不相等,这反映了不同的算法在计算效率方面是有差异的。
其中,SVD作为一种分裂法,只须按部就班地算到底,只要矩阵阶数一定,任何方程的计算量都相同。
而对于其余四种迭代法来说,不同的迭代格式在达到同一精度要求时,所需要的迭代次数、存储空间以及计算时间差异很大。
从测试中可以发现,为了达到给定的精度要求,SIRT法需要迭代55717次,而其他三种迭代法都只需迭代5次。
可以看出SIRT法的效率比较低。
综合看来,SVD法、LSQR法以及两种ICCG法在解决此类问题上的表现,远较SIRT法为优。
3 病态问题求解试验
在地震层析成像问题中,问题的稳定性是由方程（1中的矩阵A的条件数决定的,条件数越大,问题的稳定性越差。
当A的条件数非常大时,则称问题
51
2010年第10期 内蒙古石油化工
X收稿日期:
2010-03-10
是病态的。
算法的稳定性是指求解方程（1时方法的精确程度,精度越高,则算法的稳定性越好。
如此在病态问题的条件下,检验算法的稳定性。
在病态问题求解试验中,使用Hilbert矩阵来代替方程（1中的矩阵A:
Hn=[hij]∈Rn×n,hij=1
i+j-1
（i,j=1,…,n
其条件数可由以下公式计算:
k2（Hn=Kmax（Hn
Kmin（Hn
≈ean,a≈3.5
当n增加时,迅速增大,病态程度迅速增加。
本试验用n=100的Hilbert矩阵构造的方程对已有的算法程序进行测试。
为了便于比较计算结果,将解向量的各分量都设为1,再算出右侧的b。
测试结果如图1
所示。
图1 五种算法所得结果对比
从测试结果可以看出,SVD、SIRT和LSQR所得解与准确解相比差异很小,都表现出良好的稳定性。
而Jacobi-ICCG法与SSOR-ICCG法所得解跳动非常剧烈,严重偏离精确解,所以这两种算法至少在矩阵问题病态程度较高时的稳定性无法保证。
4 地质模型反演
用上述方法对给定地质模型的正演数据进行反演计算,以便进一步认识这些方法的计算效率和稳定性,为处理实际资料时选用合适的方法提供依据。
图2所示的断面模型（矩形网格划分,采用单边激发单边接收的观测系统进行地震层析成像的正演模拟,得到一组走时方程,可以写成形如（1式的稀疏矩阵方程。
再用各种算法程序求解这个方程,将所得结果与所给的断面模型进行对比,便可以看出
处理结果的正确与否。
图2 矩形分块速度模型（12行×9列和单边激发单边接收观测系统射线分布情况
图2中设定灰色网格内的速度为3km/s,深色网格内的速度为5km/s,浅色网格内的速度为2km/s。
五种方法的计算结果如图3所示。
图3 五种方法反演结果对比
从图3中可以看出SVD、SIRT和LSQR法得到的结果都比较接近所给的断面模型,而Jacobi-ICCG法和SSOR-ICCG法所得结果却与正确结果相去甚远,这个结论与图1显示的病态问题试验结果可以相互印证。
但在两种迭代法中,可以看出LSQR比SIRT的效率要高。
5 结论
综上所述,通过对SVD、SIRT、Jacobi-ICCG、SSOR-ICCG和LSQR五种可用于解决矩阵方程问题的算法的试验和比较,其中SVD、SIRT和LSQR的稳定性较好,当方程病态程度较高时,也可以得到令人满意的结果。
但SIRT法的精度低于LSQR法,效率更是远逊于后者,而SVD法在精度和效率上均居中。
所以综合比较而言,LSQR法在精度、效率和稳定性三个方面达到了较好的平衡,对于地震波走时层析成像方法而言,是合适的大型稀疏矩阵解算方法。
[参考文献]
[1] 胡家赣.线性方程组的迭代解法[M].北京:
科
学出版社,1999.
[2] 吕虹编著.矩阵论[M].北京:
科学出版社,
2001.
[3] 徐树方编著.矩阵计算的理论与方法[M].北
京:
北京大学出版社,1995.
[4] Nolet,G主编.冯锐,郝锦绮译.地震层析技术
[M].北京:
地质出版社,1991.
[5] Menke,W.著.王明光,楼海译.地球物理数据
分析——离散反演理论[M].北京:
地质出版社,1988.
[6] GeneH.Golub,CharlesF.VanLoan著.袁亚
湘等译.矩阵计算[M].北京:
科学出版社,2001.
[7] 陈景良,陈向晖著.特殊矩阵[M].北京:
清华
大学出版社,2001.
[8] 孙继广著.矩阵扰动分析[M].北京:
科学出
版社,2001.
[9] WilliamH.Press等著.傅祖芸等译.C数值
算法[M].北京:
电子工业出版社,2004.
52
内蒙古石油化工 2010年第10期
升级会员 