1、地震走時層析成像方法中大型稀疏矩陣方程解算方法比較X楊海波(成都理工大學信息工程學院,四川成都610059摘要:在地震走時層析成像方法中,對大型稀疏矩陣方程解算方的穩定性和計算效率的研究是一個重要問題。本文列舉了五種用于解決此類問題的常用算法:SVD 、SIRT 、Jacobi-ICCG 、SSOR-ICCD 和LSQR,并在求解一般問題和病態問題的試驗中,分析各方法的誤差、計算效率和穩定性;進而使用各算法對同一地質模型進行速度成像,并對結果進行比較,給出針對各種方法的客觀評價。關鍵詞:地震走時層析成像;穩定性;計算效率;誤差分析中圖分類號:P 631.4+43文獻標識碼:A 文章編號:100

2、67981(2010100051021問題的提出地震走時層析成像是根據地震波走時觀測數據求解地球內部結構的方法。在這項研究中,根據對問題的地質了解假定一個初始模型速度模型,用射線追蹤方法得到射線的理論走時,解決模型正問題,再借助于模型參數擾動,構成初始模型的理論值和觀測值之間的線性方程,反演速度結構,修改初始模型重復上述過程,直到獲得滿意的結果。一般而言,地震層析成像反演問題可以歸結為解不適定方程組:Ax=B (1通過構造正規方程:A T Ax =A Tb(2解出超定方程的最小二乘解,就可以得到這個斷面中各個單元的慢度,從而推測出斷面中不同物質的分布狀況。目前,解決這個問題主要有三類方法:&#

3、185;SVD(奇異值分解法。ºSIRT (聯立迭代圖象重建法。»投影方法。其中又有CG (共軛梯度法和LSQR (最小二乘共軛梯度法。在CG 法中,廣泛運用于實際的是PCG(預測共軛梯度法中的幾種算法,本題選用ICCG (不完全Cholisky 分解共軛梯度。本文將用上述方法求解一般問題和病態問題,分析各方法的誤差、計算效率和穩定性;進而使用它們對同一地質模型進行速度成像,并對結果進行比較,給出針對各種方法的客觀評價。2一般問題求解試驗先用規模較小的矩陣方程對各算法程序進行試運行,其目的是檢驗算法在理論上的正確性。為此構造一個系數矩陣A 的階數n=5的小型矩陣方程Ax =

4、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.0×10E-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 n

9、 =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 ,深色網格內的

11、速度為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 的穩定性較好,當方程病態程

12、度較高時,也可以得到令人滿意的結果。但SIRT 法的精度低于LSQR 法,效率更是遠遜于后者,而SVD 法在精度和效率上均居中。所以綜合比較而言,LSQR 法在精度、效率和穩定性三個方面達到了較好的平衡,對于地震波走時層析成像方法而言,是合適的大型稀疏矩陣解算方法。參考文獻1胡家贛.線性方程組的迭代解法M.北京:科學出版社,1999.2呂虹編著.矩陣論M.北京:科學出版社,2001.3徐樹方編著.矩陣計算的理論與方法M.北京:北京大學出版社,1995.4Nolet ,G 主編.馮銳,郝錦綺譯.地震層析技術M.北京:地質出版社,1991.5Menke ,W .著.王明光,樓海譯.地球物理數據分析離散反演理論M.北京:地質出版社,1988.6Ge