基于精細(xì)積分思想的反演簡單迭代算法

劉李楠, 殷興輝, 趙曉萌

(1.海軍蚌埠士官學(xué)校四系,安徽 蚌埠　233012; 2.河海大學(xué)計(jì)算機(jī)與信息學(xué)院,江蘇 南京　210098;3.安徽科技學(xué)院數(shù)理與信息工程學(xué)院,安徽 滁州　233100)

?

基于精細(xì)積分思想的反演簡單迭代算法

劉李楠1,2, 殷興輝2, 趙曉萌3

(1.海軍蚌埠士官學(xué)校四系,安徽 蚌埠233012; 2.河海大學(xué)計(jì)算機(jī)與信息學(xué)院,江蘇 南京210098;3.安徽科技學(xué)院數(shù)理與信息工程學(xué)院,安徽 滁州233100)

摘要：為了提高走時層析成像中反演算法的性能,采用一種基于精細(xì)積分的簡單迭代算法。反演計(jì)算歸結(jié)為一個簡單的迭代求解過程,對過程中出現(xiàn)的逆矩陣求解利用精細(xì)積分思想,確保迭代收斂且能收斂到方程的真解,同時具備較高的迭代速度。檢測板模型恢復(fù)測試以及實(shí)際資料反演結(jié)果表明:該方法計(jì)算過程簡單,在迭代次數(shù)較少時即能得到分辨率較高的波速剖面圖;與目前常用的一些方法相比,本文方法在反演圖像分辨率和迭代效率上都具有一定的優(yōu)勢。

關(guān)鍵詞：層析成像;反演算法;簡單迭代算法;奇異值分解法;共軛梯度法;精細(xì)積分;檢測板

人類對地球內(nèi)部物理性質(zhì)(如速度、密度、電導(dǎo)率等)以及礦產(chǎn)資源分布的了解,大多來自于可直接觀測到的地表地質(zhì)和地球物理資料的反演和解釋,由于地球物理的復(fù)雜性和不確定性,使地球物理反問題中存在許多不適定問題及病態(tài)方程組,絕大多數(shù)的觀測數(shù)據(jù)與模型參數(shù)之間都不滿足線性關(guān)系[1-5]。目前非線性反演方法主要有遺傳算法、粒子群算法及神經(jīng)網(wǎng)絡(luò)法等,這些方法的計(jì)算過程都要花費(fèi)大量的時間,同時容易受到隨機(jī)干擾的影響,對稍微復(fù)雜的介質(zhì)分布重建效果較差,所以這些方法現(xiàn)在只能用于解決小型反演問題,尚處于發(fā)展研究階段。

在一定的近似條件下,可將反演過程簡化或近似簡化為線性關(guān)系,把反演算法歸結(jié)于大型病態(tài)線性方程組的求解。有效的線性反演方法主要有:(a)奇異值分解法(SVD法)。該方法是精度較高的算法之一[6-8],缺點(diǎn)是計(jì)算速度太慢。(b)聯(lián)合代數(shù)重建算法(SIRT法)。該方法簡單、高效,但精度不高[9]。(c)共軛梯度法(LSCG法)[10]。該方法占用內(nèi)存少,但收斂不規(guī)則,與其他線性迭代方法相比,顯示出更快的收斂性及更好的可接受結(jié)果,是目前常用的反演方法?？傮w來說,利用LSCG這類迭代算法[11],解序列是最小二乘的意義上收斂的,當(dāng)連續(xù)2次迭代得到的解間距離(某種范數(shù))小于某一預(yù)先設(shè)定的值時停止迭代,顯然這樣得到的解序列難以保證收斂到問題的真解。簡單迭代算法(SI法)[12]可以解決這個問題,但其中逆矩陣的求解直接影響到計(jì)算速度以及精度。本文在一種形式上極為簡單的迭代算法基礎(chǔ)上,利用精細(xì)積分思想[13-15]求解近似逆,同時從理論上證明解序列收斂且收斂到問題的真解,提出一種基于精細(xì)積分的簡單迭代反演方法(PSI法),并通過理論模型和實(shí)際資料反演驗(yàn)證該方法的有效性。

1方 法 原 理

對于線性方程AX=B，利用如下的迭代求解過程:

(A+μI)Xk+1=B+μXk

(1)

上述迭代求解過程可以保證迭代的收斂性,且能收斂到方程的真解。但由于在地球物理反問題中,A多半是接近病態(tài)的,因此對矩陣A+μI逆的求解至關(guān)重要,將直接影響到迭代算法的精度以及計(jì)算速度。常規(guī)的高斯等方法求逆矩陣效率不高,下面利用精細(xì)積分思想求A+μI的逆矩陣。

在方程AX=B左右兩邊左乘AT,得

ATAX=ATB

(2)

H(2t)=[I+exp(Qt)]H(t)

(3)

式(3)提供了求矩陣F-1的遞推算法,接下來利用精細(xì)積分思想,取一個小步長Δt，Hm=H(2mΔt)(m=0,1,...),則

(4)

式中:m——遞推矩陣的序數(shù);j——遞推序數(shù)。

對exp(Q2iΔt)用精細(xì)積分法計(jì)算,首先對exp(QΔt)進(jìn)行Taylor展開,由于Δt是極小量,可以只保留前幾項(xiàng),得

(5)

exp(Q2iΔt)=I+Ti(i=1,2,...)

(6)

由式(3)和式(6)得

Hj+1=[I+exp(Q2iΔt)]H0(j=1,2,...)

(7)

當(dāng)j足夠大時,由式(7)即可得到精確的F-1。迭代的過程相當(dāng)于不斷倍增H(t)的積分區(qū)域,以2k的指數(shù)收斂速度收斂,因此逆矩陣求解計(jì)算效率很高。

總結(jié)上述分析,歸納出利用精細(xì)積分思想的簡單迭代反演算法(PSI)步驟如下:

步驟1用規(guī)則網(wǎng)格分割成像區(qū)域,采用平均速度法獲得初始速度模型,利用LTI法正演計(jì)算得出矩陣A,建立方程組。

步驟2選定參數(shù)μ及小步長Δt,通過式(7)計(jì)算等價方程組中ATA+μI的逆矩陣。

步驟3利用式(1)迭代計(jì)算出X,即分割后各網(wǎng)格內(nèi)慢度向量。

步驟4以步驟3中迭代得出慢度值作為新的速度模型,重復(fù)步驟1~3,直至迭代結(jié)果滿足需求。

2數(shù) 值 結(jié) 果

2.1　理論模型驗(yàn)證

定義相對誤差為

(8)

相對誤差E越小,表示所得到的估計(jì)解越接近于方程的真實(shí)解。

待測模型截面寬和深度分別為25 m和80 m,分別給出一個高速異常區(qū)和低速異常區(qū),速度分別為1 400 m/s和600 m/s,其他區(qū)域速度為1 000 m/s,位置尺寸如圖1(a)所示。采用“井間”觀測方式,地表左側(cè)井位置坐標(biāo)為(0,0),

表1　幾種方法計(jì)算效率比較Table 1　Comparison of computational efficiency of several methods

炮源置于左側(cè)井中,總炮點(diǎn)數(shù)為18個,第一個炮點(diǎn)位于深度4 m處,炮點(diǎn)均勻分布、間距為4 m,探測即檢波點(diǎn)位于右側(cè)井,檢波點(diǎn)數(shù)目及深度分布同炮點(diǎn)。

利用規(guī)則網(wǎng)格將待測區(qū)域離散化,網(wǎng)格尺度選為2.5 m×2 m,將模型區(qū)域劃分為400個小單元,假設(shè)每個小單元格內(nèi)慢度分布均勻。首先采用LTI法進(jìn)行正演走時計(jì)算,然后分別利用TSVD法、LSCG法及本文方法(PSI法)等進(jìn)行反演計(jì)算,μ值選取0.005,各種算法到達(dá)指定精度的計(jì)算效率如表1所示。

相對誤差越小越好,但并不表示相對誤差小的成像效果就一定好,在對算法進(jìn)行評價時,除了相對誤差、計(jì)算速度等定量分析,很重要的方面是直接通過反演圖像進(jìn)行定性分析。未加入噪聲,以及信噪比為13,分別利用LSCG及PSI兩種算法得出反演圖像如圖1所示(v為射線速度，以下同)。

圖1　理論模型及反演結(jié)果剖面圖Fig.1　Theoretical model and inversion results of velocity profile

可以看出,在沒有添加噪聲的情況下,2種方法反演成像都能很好地反映異常區(qū)域位置及尺寸,說明了這2種方法的有效性。進(jìn)一步仔細(xì)觀察發(fā)現(xiàn):本文PSI法描繪出的異常區(qū)域位置分布及尺寸更加符合原模型,同時其余區(qū)域速度分布更加均勻。在添加噪聲后,2種算法反演成像中異常區(qū)域都能被反映出來,高速區(qū)、低速區(qū)均明顯可見,但LSCG法成像中出現(xiàn)了一些慢度值與異常區(qū)很接近的偽跡,而利用本文PSI算法反演得到的圖慢度值與該區(qū)域理論值較接近,出現(xiàn)的少量高速干擾不會影響到異常區(qū)域的判別,基本真實(shí)地反映出模型分布情況。

進(jìn)一步分析驗(yàn)證方法的有效性,建立檢測板速度模型,觀測系統(tǒng)同上。在深度方向40 m上下分別設(shè)置不同尺寸的速度異常區(qū),低速異常區(qū)速度值為1 800 m/s,高速異常區(qū)速度值為2 200 m/s,具體分布見圖2(a)。

圖2　二維模型及其反演結(jié)果剖面圖Fig.2　Two-dimensional checkerboard model and inversion results of velocity profile

從圖2(b)、2(c)可以看出,LSCG反演圖中各區(qū)的邊緣分布不清晰,出現(xiàn)很多交錯的分布,導(dǎo)致很難正確分辨出不同的區(qū)域介質(zhì);PSI法反演圖像雖然也有少量層間交錯,但總體上邊緣分布清晰,迭代次數(shù)很少時即能基本如實(shí)描繪出不同速度區(qū)域分布。

2.2　實(shí)際數(shù)據(jù)計(jì)算

根據(jù)《河北省撫寧縣石門寨鎮(zhèn)房莊村采空區(qū)井間地震CT成像報(bào)告》[16],數(shù)據(jù)采集采用井間方式進(jìn)行,使用SWS-A多功能面波儀記錄波形,震源采用爆炸震源,電雷管激發(fā),激發(fā)點(diǎn)距3 m,接收點(diǎn)距2 m,采樣間隔0.25 ms,記錄長度256 ms。選取其中2個孔的測量數(shù)據(jù),采用平均速度法獲得初始速度模型(圖3(a)),然后利用PSI法迭代計(jì)算30次,得到波速剖面(圖3(b))。

圖3　實(shí)際數(shù)據(jù)初始速度模型及反演結(jié)果剖面圖Fig.3　Initial velocity model and real data inversion results of velocity profile

右側(cè)孔在約40 m處堵孔蓄水,40 m以下未做掃描穿透,無法取得該段相應(yīng)的測試資料,因此波速剖面圖下部存在非成像區(qū)。成像剖面速度均為2.00 ～5.15 km/s之間,淺層地表(10 m以上)波速較低(2.0～2.9 km/s之間),通常為堆積物與礫卵石層,下部為完整的砂礫巖層,波速值分布在3.35～5.16 km/s之間?；鶐r部分局部低速帶為破碎巖體段或巖溶、裂縫發(fā)育區(qū),其最大影響深度在40 m左右,40 m以下未見巖溶、破碎帶發(fā)育。以上結(jié)果表明,利用PSI法成像結(jié)果是有效的,其波速圖像直觀地反映了基巖的完整程度,與風(fēng)化溶蝕以及破碎帶分布等有著良好的對應(yīng)關(guān)系。

3結(jié)語

本文提出了一種新的地球物理線性反演算法——基于精細(xì)積分思想的簡單迭代算法,并給出了算法的具體實(shí)現(xiàn)步驟。該算法既可以保證迭代的收斂性、確保能收斂到方程的真解,同時又避免了病態(tài)或奇異矩陣求逆帶來的精度和計(jì)算速度下降,具有不依賴初值、收斂速度快等優(yōu)點(diǎn)。通過二維檢測板恢復(fù)測試比較以及實(shí)際資料反演,驗(yàn)證了PSI算法的有效性與實(shí)用性。

致謝:本文的數(shù)據(jù)由中國科學(xué)院地質(zhì)與地球物理研究所的趙連鋒老師提供,中國石油大學(xué)的黃光南博士給予了相關(guān)指導(dǎo)和建議,在此感謝他們的熱情幫助。

參考文獻(xiàn)：

[1] BENZ H M, CHOUST B A, DAWSON P B, et al. Three dimensional P and S wave velocity structure of Redoubt Volcano[J]. Journal of Geophysical Research, 1996, 101(B4): 8111-8128.

[2] HUANG Guangnan, LIU Yang. Variable damping constraint tomography and its application in VSP Data[J]. Applied Geophysics, 2012，9(2):177-185.

[3] NEMETH T, NORMARK E, QING F H. Dynamic smoothing in cross well travel time tomography[J]. Geophysics, 1997,62(1):168-176.

[4] MAURER H, GREEN A G, HOLLIGER K. Full-waveform inversion of cross-hole radar data based on 2-D finite-difference time-domain solutions of Maxwell’s equations[J]. IEEE Transactions on Geo-Science and Remote Sensing, 2007,72(5):53-64.

[5] PIDLISECKY A, HABER E, KNIGHT R. RESINVM3D: a 3D resistivity inversion package[J]. Geophysics, 2007, 72(2): 1-10.

[6] PAIGE C C, SAUNDER M A. Sparse linear equations and least squares problems[J]. ACM Transactions on Mathematical Software, 1982, 8(2):195-209.

[7] LINDE N A, TRYGGYASON J E, HUBBARD S S, et al. Joint inversion of cross hole radar and seismic travel times acquired at the south oyster bacterial transport site[J]. Geophysics, 2008,73(4):39-50.

[8] 顧勇為,歸慶明,邊少鋒. 地球物理反問題中兩種正則化方法的比較[J]. 武漢大學(xué)學(xué)報(bào)：信息科學(xué)版,2005,30(3):238-241. (GU Yongwei, GUI Qingming, BIAN Shaofeng. Comparison of two kinds of regularization methods in geophysical inverse problem[J]. Journal of Wuhan University：Information Science Edition, 2005,30(3):238-241.(in Chinese))

[9] 曹俊興,聶在平. 基于物理模型的成像算法評價與SASART成像算法[J]. 成都理工學(xué)院學(xué)報(bào), 1998,25(4):473-479.(CAO Junxing, NIE Zaiping. Evaluation of imaging algorithm and SASART imaging algorithm based on the physical model[J]. Journal of Chengdu Institute of Science and Technology, 1998, 25(4):473-479. (in Chinese))

[10] SLEIPEN G L, FOKKEMA D R. Bicgstab(L) for linear equations involving unsymmetric matrices with complex spectrum[J].Electronic Trasactions on Numerical Analysis,1993, 1(11):11-32.

[11] RAWLINSONA N, PZAGAYA S, FISHWICHKB S. Seismic tomography: a window into deep Earth[J]. Physics of the Earth and Planetary Interiors, 2010,178(3):101-135.

[12] 毛先進(jìn),楊玲英. 病態(tài)線性方程組的簡單迭代解法[J]. 物探化探計(jì)算技術(shù),1999, 21(1):14-18. (MAO Xianjin, YANG Lingying. A simple iterative method for solving ill conditioned linear system of equations[J]. Computing Techniques for Geophysical and Geochemical Exploration,1999, 21(1):14-18. (in Chinese))

[13] 富明慧,張文志. 病態(tài)線性方程的精細(xì)積分解法[J]. 計(jì)算力學(xué)學(xué)報(bào),2011, 28(4):530-534. (FU Minghui, ZHANG Wenzhi. Precise integration method for solving ill conditioned linear equations[J]. Chinese Journal of Computational Mechanics, 2011,28(4):530-534. (in Chinese))

[14] 張文志,黃培彥. 病態(tài)代數(shù)系統(tǒng)求解的精細(xì)迭代方法[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2013, 34(7):736-741. (ZHANG Wenzhi, HUANG Peiyan. Precise iterative method for solving ill conditioned algebraic system[J]. Applied Mathematics and Mechanics, 2013, 34(7):736-741. (in Chinese))

[15] 楊綠峰,洪斌,高欽,等. 混凝土結(jié)構(gòu)中氯離子擴(kuò)散分析的精細(xì)積分法[J]. 水利水電科技進(jìn)展,2012, 32(2): 32-36. (YANG Lyufeng, HONG Bin, GAO Qin, et al. Precise integration method for analysis of chloride diffusion in mix soil structure[J]. Progress of Science and Technology of Water Resources, 2012, 32 (2): 32-36. (in Chinese))

[16] 李洪濤,孫黨生. 河北省撫寧縣石門寨鎮(zhèn)房莊村采空區(qū)井間地震CT成像報(bào)告[R].保定：地礦保定工程勘察院, 2000.

Simple iterative inversion algorithm based on precise integration method

LIU Linan1,2, YIN Xinghui2, ZHAO Xiaomeng3

(1.TheFourthDepartment,BengbuNavalPettyOfficerAcademy,Bengbu233012,China;

2.CollegeofComputerandInformation,HohaiUniversity,Nanjing210098,China;

3.CollegeofMathematics,PhysicsandInformationEngineering,AnhuiScienceand

TechnologyUniversity,Chuzhou233100,China)

Abstract:In order to improve the performance of the inversion algorithm for travel time tomography, a simple iterative algorithm based on the precise integration method was used. The improved algorithm could be considered a simple iterative process. The precise integration method was used to solve the inverse matrix to ensure convergence and obtain the true solution to the equation at a higher speed of iteration. The checkerboard model recovery test and real data inversion results show that the calculation process of the improved algorithm is simple, and the algorithm can obtain the velocity profile with higher resolution and fewer iterations. Compared with some commonly used methods, this algorithm has higher image resolution and greater iteration efficiency.

Key words:tomography; inversion algorithm; simple iterative algorithm; singular value decomposition; conjugate gradient method; precise integration; checkerboard

中圖分類號：P631

文獻(xiàn)標(biāo)志碼：A

文章編號：1000-1980(2015)06-0594-05

作者簡介：劉李楠(1983—),男,安徽六安人,博士研究生,主要從事電磁波及聲波層析成像研究。E-mail: linan10241@sina.com

基金項(xiàng)目：安徽省級高校自然科學(xué)基金(KJ2012Z 056);安徽科技學(xué)院青年基金(ZRC2014442)

收稿日期：2014-12-10

DOI：10.3876/j.issn.1000-1980.2015.06.015