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

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

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

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

波速層析成像的概念可以由下面這個簡單的數(shù)學公式來表示:

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

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

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

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

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

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

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

1.2 Dijkstra最短路徑射線追蹤法

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

1.3 建立JACOBI矩陣

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

1.4 反演計算

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

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

2 假定物理模型

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

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

模型1(圖2)設置水平低速異常帶，帶寬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)設置4個速度異常塊，2個高速異常塊，2個低速異常塊，尺寸均為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

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

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

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

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

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

3 結論

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

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

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

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

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

［3］趙明階，許錫賓.水工混凝土結構缺陷成像診斷的試驗研究［J］.水利學報，2007，38(2):198 －204.

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

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

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

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

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

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