基于最小走時(shí)射線的LSQR 成像反演方法研究

趙火焱，趙明階，黃衛(wèi)東，李 建，李 勇

(1.重慶交通大學(xué)河海學(xué)院，重慶400074;2.重慶市交通工程監(jiān)理咨詢(xún)有限責(zé)任公司，重慶400042)

層析成像作為一種非破壞性的檢測(cè)方法，概念是利用外部探測(cè)能量觀察待測(cè)物的反應(yīng)，藉由反應(yīng)結(jié)果進(jìn)而對(duì)待測(cè)物剖面進(jìn)行成像的技術(shù)。其最早應(yīng)用于天文顯像(Bracewell，1956)、醫(yī)學(xué)診斷與顯微技術(shù)(DeRosier and Klug，1968)等，并逐漸應(yīng)用于地球物理學(xué)和機(jī)械工業(yè)中。如今隨著計(jì)算機(jī)技術(shù)和成像理論的發(fā)展，成像技術(shù)在土木工程無(wú)損檢測(cè)中也日益發(fā)揮重要作用［1－3］。采用彈性波等激發(fā)源，根據(jù)實(shí)測(cè)波初至旅行時(shí)間及波的行進(jìn)路徑來(lái)反演待測(cè)體內(nèi)速度分布的方法稱(chēng)為波速層析成像。

波速層析成像的概念可以由下面這個(gè)簡(jiǎn)單的數(shù)學(xué)公式來(lái)表示:

為確定波的行進(jìn)路線ray，目前主要采用射線追蹤技術(shù)，其理論基礎(chǔ)是在高頻近似條件下，波場(chǎng)的主能量沿射線軌跡傳播?，F(xiàn)有射線追蹤方法有初至波旅行時(shí)法、線性旅行時(shí)插值法和最短路徑法等。其中最短路徑射線追蹤算法最早由Nakanishi et al提出，Moser［4］對(duì)此方法做了全面研究。該方法基于圖論和費(fèi)馬原理，其靈活而穩(wěn)定，且適用于任意復(fù)雜介質(zhì)模型，在彈性波正反演領(lǐng)域應(yīng)用廣泛［5－6］。

筆者根據(jù)最短路徑射線追蹤原理，并在離散網(wǎng)格單元過(guò)程中做出改進(jìn)，采用FORTRAN語(yǔ)言和VB語(yǔ)言實(shí)現(xiàn)聯(lián)合編程計(jì)算，然后通過(guò)建立數(shù)字模型來(lái)驗(yàn)證該方法的可行性。

1 波速層析成像系統(tǒng)的軟件編制

1.1 建立速度網(wǎng)格和射線網(wǎng)格

為實(shí)現(xiàn)編程計(jì)算，筆者將待測(cè)體離散化為N個(gè)二維矩形節(jié)點(diǎn)網(wǎng)格作為速度網(wǎng)格，并假定網(wǎng)格單元內(nèi)速度v為均勻分布。第i個(gè)單元的初始速度值由旅行時(shí)觀測(cè)值假定:該方法首先假定模型為均質(zhì)，則射線傳播路徑為直線，由第j條射線的長(zhǎng)度Dj為與這條射線對(duì)應(yīng)的觀測(cè)值的比值得到射線上的平均速度值，令通過(guò)某網(wǎng)格單元的所有射線的速度均值為該單元的初始速度值。該方法比假設(shè)初始速度值為均值的方法具有一定的優(yōu)越性［7］。

將速度網(wǎng)格各邊內(nèi)插若干節(jié)點(diǎn)并與原節(jié)點(diǎn)網(wǎng)格共同形成射線網(wǎng)格。本文用到的射線網(wǎng)格模型是在Moser使用的網(wǎng)絡(luò)模型的基礎(chǔ)上每單元增加了4個(gè)角點(diǎn)，即每條邊增加2個(gè)節(jié)點(diǎn)，這在一定程度上提高了射線精度。根據(jù)惠更斯原理，每個(gè)射線節(jié)點(diǎn)既作為接收點(diǎn)又作為新的震源向周?chē)?jié)點(diǎn)傳播能量，波在各節(jié)點(diǎn)之間的旅行時(shí)由下式得出:

其中:tij、dij為節(jié)點(diǎn)間走時(shí)及間距;Sij為節(jié)點(diǎn)公共單元慢度。

1.2 Dijkstra最短路徑射線追蹤法

求解最短路徑問(wèn)題多采用F.W.Dijkstra［8］提出Dijkstra算法。在求從網(wǎng)絡(luò)中的某一節(jié)點(diǎn)(源點(diǎn))到其余各節(jié)點(diǎn)的最短路徑時(shí)，該算法將網(wǎng)絡(luò)中的節(jié)點(diǎn)分成3部分:未標(biāo)記節(jié)點(diǎn)、臨時(shí)標(biāo)記節(jié)點(diǎn)和最短路徑節(jié)點(diǎn)。算法開(kāi)始時(shí)源點(diǎn)初始化為最短路徑節(jié)點(diǎn)，其余為未標(biāo)記節(jié)點(diǎn)，算法執(zhí)行過(guò)程中，每次從最短路徑節(jié)點(diǎn)往相鄰節(jié)點(diǎn)擴(kuò)展，非最短路徑節(jié)點(diǎn)的相鄰節(jié)點(diǎn)修改為臨時(shí)標(biāo)記節(jié)點(diǎn)，判斷權(quán)值是否更新后，在所有臨時(shí)標(biāo)記節(jié)點(diǎn)中提取權(quán)值最小的節(jié)點(diǎn)，修改為最短路徑節(jié)點(diǎn)后作為下一次的擴(kuò)展源，再重復(fù)前面的步驟，當(dāng)所有節(jié)點(diǎn)都做過(guò)擴(kuò)展源后算法結(jié)束。具體算法描述如下:設(shè)在一非負(fù)權(quán)簡(jiǎn)單連通無(wú)向圖G=［V，E，W ］(V:頂點(diǎn)集;E:邊集;W:邊權(quán)值)中，D 為圖G的鄰接矩陣，求源點(diǎn)P0到其余所有節(jié)點(diǎn)P(i)的最短路徑長(zhǎng)度。①將V分為未標(biāo)記節(jié)點(diǎn)子集N、臨時(shí)最短路徑節(jié)點(diǎn)子集T和最短路徑節(jié)點(diǎn)子集S，每個(gè)節(jié)點(diǎn)上的路徑權(quán)值為 D(i)，初始化:S=P0，T=φ，N=V－C，D(0)=O，D(i)=∞;②更新:將新加入 S集合的節(jié)點(diǎn)P(s)作為擴(kuò)展源，計(jì)算從擴(kuò)展源到相鄰節(jié)點(diǎn)的路徑值。若該值比原值小，則替換原值，否則保持原值不變，將相鄰節(jié)點(diǎn)之中的未標(biāo)記節(jié)點(diǎn)歸為臨時(shí)標(biāo)記節(jié)點(diǎn)，即T=T∪ P(i)，N=N－ P(i);③選擇:在T中選擇具有最小路徑值D(s)的節(jié)點(diǎn)P(s)，歸入集合S中，即 S=S∪P(s)，T=T－P(s);④迭代判斷:若T=φ算法結(jié)束，否則轉(zhuǎn)②。

1.3 建立JACOBI矩陣

JACOBI矩陣A(M，N)是反映速度節(jié)點(diǎn)與各條射線關(guān)系的矩陣，M為射線根數(shù)，N為速度網(wǎng)格單元個(gè)數(shù)。當(dāng)最短射線路徑確定后，A(i，j)即為第i條射線在第j單元中的長(zhǎng)度。由于每條射線只通過(guò)少數(shù)幾個(gè)單元，只對(duì)這幾個(gè)單元的節(jié)點(diǎn)產(chǎn)生影響，故JACOBI矩陣是大型稀疏矩陣，并且通常是病態(tài)的。

1.4 反演計(jì)算

網(wǎng)格離散化后，第j條射線旅行時(shí)理論值Tt

迭代方法有代數(shù)重建法(ART)、聯(lián)合迭代重建法(SIRT)以及最小二乘正交分解法(LSQR)等幾種。LSQR算法［9］計(jì)算速度快、易收斂，故采用此種方法。

2 假定物理模型

筆者設(shè)計(jì)了2個(gè)帶有異常速度區(qū)域的假定物理模型。模型寬6 m，深8 m，離散網(wǎng)格間距0.5 m，每邊內(nèi)插9個(gè)節(jié)點(diǎn)，發(fā)射點(diǎn)22個(gè)，接收點(diǎn)24個(gè)，炮點(diǎn)間距均為0.5 m，左側(cè)發(fā)射時(shí)頂部接收，頂部發(fā)射時(shí)左右兩側(cè)接收，射線布置圖如圖1。

圖1 射線布置圖Fig.1 Ray disposal figure

模型1(圖2)設(shè)置水平低速異常帶，帶寬2 m，位置在頂面以下2 m，背景速度4 000 m/s，異常帶速度為3 500 m/s。經(jīng)7次迭代后射線圖和色塊圖及波速三維圖分別為圖3～圖5。

圖2 模型1Fig.2 Model No.1

圖3 經(jīng)7次迭代射線圖Fig.3 Ray after 7 times inversion

圖4 經(jīng)7次迭代后色塊圖Fig.4 Result after 7 times inversion

圖5 經(jīng)7次迭代后波速3D圖Fig.5 Velocity 3D figure after 7 times inversion

模型2(圖6)設(shè)置4個(gè)速度異常塊，2個(gè)高速異常塊，2個(gè)低速異常塊，尺寸均為1 m×1 m。背景速度4 000 m/s，高低速異常塊速度分別為4 500，3 500 m/s。經(jīng)6次迭代后的射線圖和色塊圖及波速三維圖分別為圖7～圖9。

圖6 模型2Fig.6 Model No.2

圖7 經(jīng)6次迭代后射線圖Fig.7 Ray after 6 times inversion

圖8 經(jīng)6次迭代后色塊圖Fig.8 Result after 6 times inversion

圖9 經(jīng)6次迭代后波速3D圖Fig.9 Velocity 3D figure after 6 times inversion

通過(guò)觀察圖4及圖8可見(jiàn)，成像結(jié)果能較好的反映速度異常區(qū)域的位置，但在形狀和尺寸上存在失真，究其原因有以下幾點(diǎn):①本文模型中，射線是水平和傾斜的，缺少垂直方向或接近垂直方向上的射線，因此成像橫向分辨率較垂向分辨率要高;②由于通過(guò)底部區(qū)域的射線數(shù)較少，導(dǎo)致這一區(qū)域的失真情況較明顯;③由于在邊界處數(shù)據(jù)的覆蓋程度不足引起邊界圖像失真并影響反演的中間區(qū)域，使遠(yuǎn)離邊界很大一段距離內(nèi)的圖象都發(fā)生畸變。

模型節(jié)點(diǎn)速度真值與反演值之間相關(guān)系數(shù)為:

模型1反演結(jié)果的平均相對(duì)誤差為2.34%，相關(guān)系數(shù)為0.84，低速異常帶相對(duì)誤差為－4.4% ～12.5%;模型2平均相對(duì)誤差為2.59%，相關(guān)系數(shù)為0.85，低速區(qū)節(jié)點(diǎn)相對(duì)誤差為 －1.73% ～7.44%，高速區(qū)相對(duì)誤差為－6.24% ～4.06%。節(jié)點(diǎn)相對(duì)誤差如圖10、圖11。結(jié)果表明，異常區(qū)偏差較大，且結(jié)果值大都偏近于速度背景值;平均相對(duì)誤差較小，相關(guān)系數(shù)值較高;成像結(jié)果能較好地反映模型值。

圖10 模型1節(jié)點(diǎn)相對(duì)誤差Fig.10 Node relative difference of model NO.1

圖11 模型2節(jié)點(diǎn)相對(duì)誤差Fig.11 Node relative difference of model NO.2

3 結(jié)論

對(duì)基于最小走時(shí)射線追蹤方法和LSQR反演方法的波速層析成像從方法原理、假定物理模型等方面進(jìn)行了研究，得出了以下2點(diǎn)結(jié)論:

1)基于最小走時(shí)射線追蹤的LSQR波速層析成像方法可靠穩(wěn)定，收斂性好，效率高。

2)由于射線分布不均及射線數(shù)據(jù)覆蓋程度不足等原因引起圖像局部失真。以上情況可通過(guò)完善射線布置情況和加密網(wǎng)格等方法得以改善。

［1］楊文采，杜劍淵.層析成像新算法及其在工程檢測(cè)上的應(yīng)用［J］.地球物理學(xué)報(bào)，1994，37(2):239 －244.

［2］王五平，宋人心，傅翔，等.用超聲波CT探測(cè)混凝土內(nèi)部缺陷［J］.水利水運(yùn)工程學(xué)報(bào)，2003，56(5):56 －60.

［3］趙明階，許錫賓.水工混凝土結(jié)構(gòu)缺陷成像診斷的試驗(yàn)研究［J］.水利學(xué)報(bào)，2007，38(2):198 －204.

［4］Moser T J.Shortest path calculation of seismic rays［J］.Geophysics，1991，56(1):59 －67.

［5］張建中，陳世軍，余大祥.最短路徑射線追蹤方法及其改進(jìn)［J］.地球物理學(xué)進(jìn)展，2003(1):146－150.

［6］周兵，趙明階.最小走時(shí)射線追蹤層析方法［J］.物化探計(jì)算技術(shù)，1992(2):124－130.

［7］段心標(biāo)，金維浚.井間地震層析成像初始速度模型［J］.地球物理學(xué)進(jìn)展，2007(6):1831－1835.

［8］王志和，凌云.Dijkstra最短路徑算法的優(yōu)化及其實(shí)現(xiàn)［J］.微計(jì)算機(jī)信息，2007(33):275－277.

［9］常旭，盧孟夏，劉伊克.地震層析成像反演中3種廣義解的誤差分析與評(píng)價(jià)［J］.地球物理學(xué)報(bào)，1999(5):695－702.