地震走时层析成像方法中大型稀疏矩阵方程解算方法比较.docx

1、地震走时层析成像方法中大型稀疏矩阵方程解算方法比较地震走时层析成像方法中大型稀疏矩阵方程解算方法比较X杨海波(成都理工大学信息工程学院,四川成都610059摘要:在地震走时层析成像方法中,对大型稀疏矩阵方程解算方的稳定性和计算效率的研究是一个重要问题。本文列举了五种用于解决此类问题的常用算法:SVD 、SIRT 、Jacobi-ICCG 、SSOR-ICCD 和LSQR,并在求解一般问题和病态问题的试验中,分析各方法的误差、计算效率和稳定性;进而使用各算法对同一地质模型进行速度成像,并对结果进行比较,给出针对各种方法的客观评价。关键词:地震走时层析成像;稳定性;计算效率;误差分析中图分类号:P

2、 631.4+43文献标识码:A 文章编号:10067981(2010100051021问题的提出地震走时层析成像是根据地震波走时观测数据求解地球内部结构的方法。在这项研究中,根据对问题的地质了解假定一个初始模型速度模型,用射线追踪方法得到射线的理论走时,解决模型正问题,再借助于模型参数扰动,构成初始模型的理论值和观测值之间的线性方程,反演速度结构,修改初始模型重复上述过程,直到获得满意的结果。一般而言,地震层析成像反演问题可以归结为解不适定方程组:Ax=B (1通过构造正规方程:A T Ax =A Tb(2解出超定方程的最小二乘解,就可以得到这个断面中各个单元的慢度,从而推测出断面中不同物质

3、的分布状况。目前,解决这个问题主要有三类方法:SVD(奇异值分解法。SIRT (联立迭代图象重建法。投影方法。其中又有CG (共轭梯度法和LSQR (最小二乘共轭梯度法。在CG 法中,广泛运用于实际的是PCG(预测共轭梯度法中的几种算法,本题选用ICCG (不完全Cholisky 分解共轭梯度。本文将用上述方法求解一般问题和病态问题,分析各方法的误差、计算效率和稳定性;进而使用它们对同一地质模型进行速度成像,并对结果进行比较,给出针对各种方法的客观评价。2一般问题求解试验先用规模较小的矩阵方程对各算法程序进行试运行,其目的是检验算法在理论上的正确性。为此构造一个系数矩阵A 的阶数n=5的小型矩

4、阵方程Ax =b(按方程(2,A 采用对称正定矩阵的形式:2.42.56.94.31.32.55.66.12.87.76.96.14.98.15.74.32.88.15.91.81.37.75.71.83.5x 1x 2x 3x 4x 5=4.65.16.52.43.3该方程的准确解为:X =(- 4.635068,0.279888,-0.670548,4.287092,0.935946T ,各个方法算出的解为x ,设d =Ax -b 定义了解出解与准确解之间的差向量。表1列出了准确解和SV D 、SIRT 、Jacobi-ICCG 、SSOR -ICCD 、LSQR 五种方法得到的解和误差,

5、迭代计算在解的残差小于1.010E-6时停止计算。表1五种方法解算结果对比准确解SVDSIRT (迭代次数:55717Jacobi -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.287092 4.287092-0.000000

6、 4.287092-0.000000 4.2870890.0000010.9359460.9359460.0000000.9359460.0000000.9359460.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.287092 4.287092-0.000000

7、 4.287092-0.0000000.9359460.9359460.0000000.9359460.000000从表1可以看出,所有的五个算法能得到非常接近准确值的解,这说明所有的算法在理论上都是正确的。但很明显地,五个算法在得到相近结果的同时,所付出的计算上的代价并不相等,这反映了不同的算法在计算效率方面是有差异的。其中,SVD 作为一种分裂法,只须按部就班地算到底,只要矩阵阶数一定,任何方程的计算量都相同。而对于其余四种迭代法来说,不同的迭代格式在达到同一精度要求时,所需要的迭代次数、存储空间以及计算时间差异很大。从测试中可以发现,为了达到给定的精度要求,SIRT 法需要迭代55717

8、次,而其他三种迭代法都只需迭代5次。可以看出SIRT 法的效率比较低。综合看来,SVD 法、LSQR 法以及两种ICCG 法在解决此类问题上的表现,远较SIRT 法为优。3病态问题求解试验在地震层析成像问题中,问题的稳定性是由方程(1中的矩阵A 的条件数决定的,条件数越大,问题的稳定性越差。当A 的条件数非常大时,则称问题512010年第10期内蒙古石油化工X 收稿日期:2010-03-10是病态的。算法的稳定性是指求解方程(1时方法的精确程度,精度越高,则算法的稳定性越好。如此在病态问题的条件下,检验算法的稳定性。在病态问题求解试验中,使用Hilbert 矩阵来代替方程(1中的矩阵A :H

9、n =h ij R n n ,h ij =1i+j-1(i,j=1,n其条件数可由以下公式计算:k 2(H n =K max (H n K min (H n e an ,a 3.5当n 增加时,迅速增大,病态程度迅速增加。本试验用n =100的Hilbert 矩阵构造的方程对已有的算法程序进行测试。为了便于比较计算结果,将解向量的各分量都设为1,再算出右侧的b 。测试结果如图1 所示。图1五种算法所得结果对比从测试结果可以看出,SVD 、SIRT 和LSQR 所得解与准确解相比差异很小,都表现出良好的稳定性。而Jacobi -ICCG 法与SSOR -ICCG 法所得解跳动非常剧烈,严重偏离精

10、确解,所以这两种算法至少在矩阵问题病态程度较高时的稳定性无法保证。4地质模型反演用上述方法对给定地质模型的正演数据进行反演计算,以便进一步认识这些方法的计算效率和稳定性,为处理实际资料时选用合适的方法提供依据。图2所示的断面模型(矩形网格划分,采用单边激发单边接收的观测系统进行地震层析成像的正演模拟,得到一组走时方程,可以写成形如(1式的稀疏矩阵方程。再用各种算法程序求解这个方程,将所得结果与所给的断面模型进行对比,便可以看出 处理结果的正确与否。图2矩形分块速度模型(12行9列和单边激发单边接收观测系统射线分布情况图2中设定灰色网格内的速度为3km /s ,深色网格内的速度为5km/s,浅色

11、网格内的速度为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 的稳定性较好,当方程病态程度较高时,也可以得到令

12、人满意的结果。但SIRT 法的精度低于LSQR 法,效率更是远逊于后者,而SVD 法在精度和效率上均居中。所以综合比较而言,LSQR 法在精度、效率和稳定性三个方面达到了较好的平衡,对于地震波走时层析成像方法而言,是合适的大型稀疏矩阵解算方法。参考文献1胡家赣.线性方程组的迭代解法M.北京:科学出版社,1999.2吕虹编著.矩阵论M.北京:科学出版社,2001.3徐树方编著.矩阵计算的理论与方法M.北京:北京大学出版社,1995.4Nolet ,G 主编.冯锐,郝锦绮译.地震层析技术M.北京:地质出版社,1991.5Menke ,W .著.王明光,楼海译.地球物理数据分析离散反演理论M.北京:地质出版社,1988.6Gene H.Golub,Charles F.Van Loan 著.袁亚湘等译.矩阵计算M .北京:科学出版社,2001.7陈景良,陈向晖著.特殊矩阵M .北京:清华大学出版社,2001.8孙继广著.矩阵扰动分析M .北京:科学出版社,2001.9William H .Press 等著.傅祖芸等译.C 数值算法M.北京:电子工业出版社,2004.52内蒙古石油化工2010年第10期