地震走時(shí)層析成像方法中大型稀疏矩陣方程解算方法比較

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

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

3、185;SVD(奇異值分解法。ºSIRT (聯(lián)立迭代圖象重建法。»投影方法。其中又有CG (共軛梯度法和LSQR (最小二乘共軛梯度法。在CG 法中,廣泛運(yùn)用于實(shí)際的是PCG(預(yù)測共軛梯度法中的幾種算法,本題選用ICCG (不完全Cholisky 分解共軛梯度。本文將用上述方法求解一般問題和病態(tài)問題,分析各方法的誤差、計(jì)算效率和穩(wěn)定性;進(jìn)而使用它們對同一地質(zhì)模型進(jìn)行速度成像,并對結(jié)果進(jìn)行比較,給出針對各種方法的客觀評價(jià)。2一般問題求解試驗(yàn)先用規(guī)模較小的矩陣方程對各算法程序進(jìn)行試運(yùn)行,其目的是檢驗(yàn)算法在理論上的正確性。為此構(gòu)造一個(gè)系數(shù)矩陣A 的階數(shù)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該方程的準(zhǔn)確解為:X =(- 4.635068,0.279888,-0.670548,4.287092,0.935946T ,各個(gè)方法算出的解為x ,設(shè)d =Ax -b 定義了解出解與準(zhǔn)確解之間的差向量。表1列出了準(zhǔn)確解和SV D 、SIRT 、Jacobi-ICCG 、SSOR -ICCD 、LSQR 五種方法得到的解和誤差,迭代計(jì)算在解的

5、殘差小于1.0×10E-6時(shí)停止計(jì)算。表1五種方法解算結(jié)果對比準(zhǔn)確解SVDSIRT (迭代次數(shù):55717Jacobi -ICCG (迭代次數(shù):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準(zhǔn)確解SSOR-ICCG (迭代次數(shù):5,松弛因子:X =1.0LSQR (迭代次數(shù):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可以看出,所有的五個(gè)算法能得到非常接近準(zhǔn)確值的解,這說明所有的算法在理論上都是正確的。但很明顯地,五個(gè)算法在得到相近結(jié)果的同時(shí),所付出的計(jì)算上的代價(jià)并不相等,這反映了不同的算法在計(jì)算效率方面是有差異的。其中,SVD 作為一種分裂法,只須按部就班地算到底,只要矩陣階數(shù)一定,任何方程的計(jì)算量都相同。而對于其余四種迭代法來說,不同的迭代格式在達(dá)到同一精度要求時(shí),所需要的迭代次數(shù)、存儲空間以及計(jì)算時(shí)間差異很大。從測試中可以發(fā)現(xiàn),為了達(dá)到給定的精度要求,SIRT 法需要迭代55717次

8、,而其他三種迭代法都只需迭代5次。可以看出SIRT 法的效率比較低。綜合看來,SVD 法、LSQR 法以及兩種ICCG 法在解決此類問題上的表現(xiàn),遠(yuǎn)較SIRT 法為優(yōu)。3病態(tài)問題求解試驗(yàn)在地震層析成像問題中,問題的穩(wěn)定性是由方程(1中的矩陣A 的條件數(shù)決定的,條件數(shù)越大,問題的穩(wěn)定性越差。當(dāng)A 的條件數(shù)非常大時(shí),則稱問題512010年第10期內(nèi)蒙古石油化工X 收稿日期:2010-03-10是病態(tài)的。算法的穩(wěn)定性是指求解方程(1時(shí)方法的精確程度,精度越高,則算法的穩(wěn)定性越好。如此在病態(tài)問題的條件下,檢驗(yàn)算法的穩(wěn)定性。在病態(tài)問題求解試驗(yàn)中,使用Hilbert 矩陣來代替方程(1中的矩陣A :H n

9、 =h ij R n ×n ,h ij =1i+j-1(i,j=1,n其條件數(shù)可由以下公式計(jì)算:k 2(H n =K max (H n K min (H n e an ,a 3.5當(dāng)n 增加時(shí),迅速增大,病態(tài)程度迅速增加。本試驗(yàn)用n =100的Hilbert 矩陣構(gòu)造的方程對已有的算法程序進(jìn)行測試。為了便于比較計(jì)算結(jié)果,將解向量的各分量都設(shè)為1,再算出右側(cè)的b 。測試結(jié)果如圖1 所示。圖1五種算法所得結(jié)果對比從測試結(jié)果可以看出,SVD 、SIRT 和LSQR 所得解與準(zhǔn)確解相比差異很小,都表現(xiàn)出良好的穩(wěn)定性。而Jacobi -ICCG 法與SSOR -ICCG 法所得解跳動非常劇烈,

10、嚴(yán)重偏離精確解,所以這兩種算法至少在矩陣問題病態(tài)程度較高時(shí)的穩(wěn)定性無法保證。4地質(zhì)模型反演用上述方法對給定地質(zhì)模型的正演數(shù)據(jù)進(jìn)行反演計(jì)算,以便進(jìn)一步認(rèn)識這些方法的計(jì)算效率和穩(wěn)定性,為處理實(shí)際資料時(shí)選用合適的方法提供依據(jù)。圖2所示的斷面模型(矩形網(wǎng)格劃分,采用單邊激發(fā)單邊接收的觀測系統(tǒng)進(jìn)行地震層析成像的正演模擬,得到一組走時(shí)方程,可以寫成形如(1式的稀疏矩陣方程。再用各種算法程序求解這個(gè)方程,將所得結(jié)果與所給的斷面模型進(jìn)行對比,便可以看出 處理結(jié)果的正確與否。圖2矩形分塊速度模型(12行×9列和單邊激發(fā)單邊接收觀測系統(tǒng)射線分布情況圖2中設(shè)定灰色網(wǎng)格內(nèi)的速度為3km /s ,深色網(wǎng)格內(nèi)的

11、速度為5km/s,淺色網(wǎng)格內(nèi)的速度為2km/s 。 五種方法的計(jì)算結(jié)果如圖3所示。圖3五種方法反演結(jié)果對比從圖3中可以看出SVD 、SIRT 和LSQR 法得到的結(jié)果都比較接近所給的斷面模型,而Jacobi -ICCG 法和SSOR-ICCG 法所得結(jié)果卻與正確結(jié)果相去甚遠(yuǎn),這個(gè)結(jié)論與圖1顯示的病態(tài)問題試驗(yàn)結(jié)果可以相互印證。但在兩種迭代法中,可以看出LSQR 比SIRT 的效率要高。5結(jié)論綜上所述,通過對SVD 、SIRT 、Jacobi -ICCG 、SSOR -ICCG 和LSQR 五種可用于解決矩陣方程問題的算法的試驗(yàn)和比較,其中SVD 、SIRT 和LSQR 的穩(wěn)定性較好,當(dāng)方程病態(tài)程

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