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

趙火焱，趙明階，黃衛東，李 建，李 勇

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

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

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

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

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

1 波速層析成像系統的軟件編制

1.1 建立速度網格和射線網格

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

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

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

1.2 Dijkstra最短路徑射線追蹤法

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

1.3 建立JACOBI矩陣

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

1.4 反演計算

網格離散化后，第j條射線旅行時理論值Tt

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

2 假定物理模型

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

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

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

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

圖3 經7次迭代射線圖Fig.3 Ray after 7 times inversion

圖4 經7次迭代后色塊圖Fig.4 Result after 7 times inversion

圖5 經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。經6次迭代后的射線圖和色塊圖及波速三維圖分別為圖7～圖9。

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

圖7 經6次迭代后射線圖Fig.7 Ray after 6 times inversion

圖8 經6次迭代后色塊圖Fig.8 Result after 6 times inversion

圖9 經6次迭代后波速3D圖Fig.9 Velocity 3D figure after 6 times inversion

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

模型節點速度真值與反演值之間相關系數為:

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

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

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

3 結論

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

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

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

［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最短路徑算法的優化及其實現［J］.微計算機信息，2007(33):275－277.

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