堤防CT成像的數(shù)值模擬研究

潘紀順 劉宇鋒 李長征

摘 要：簡述以射線理論為基礎(chǔ)的地震CT的基本原理、數(shù)學方法、圖像重構(gòu)技術(shù)、ART算法和SIRT算法，在現(xiàn)有井間地震CT的基礎(chǔ)上研究堤防CT，利用MATLAB編程對堤防進行地震斷面CT的數(shù)值模擬，通過ART和SIRT重建方法進行反演，分析得出二者在堤防CT應(yīng)用中的優(yōu)缺點及適用條件。結(jié)果表明：ART算法和SIRT算法都得出了較好的反演結(jié)果，提高了探測效率和精度，較好地驗證了堤防CT成像技術(shù)的正確性和可行性。

關(guān)鍵詞：堤防;地震CT;ART算法;SIRT算法

Abstract： The basic principle， mathematical method， image reconstruction technology， ART algorithm and SIRT algorithm of seismic CT based on ray theory were briefly introduced. Based on the existing cross well CT， the embankment CT was studied and the seismic cross-section CT numerical simulation was carried out by MATLAB programming. The advantages and disadvantages of ART and SIRT reconstruction methods in the application of embankment CT were analyzed. Finally， good inversion results were obtained， which improved the accuracy and verified the correctness and feasibility of embankment CT imaging technology.

Key words： embankment; seismic CT; ART algorithm; SIRT algorithm

1 引 言

長久以來，每到汛期，洪災(zāi)威脅對兩岸居民的安全和生活造成了嚴重的影響。抵御洪水泛濫、保護人民生活和工農(nóng)業(yè)的生產(chǎn)，堤防已經(jīng)成為一道重要的屏障。我國江河堤防質(zhì)量良莠不齊，為了有效防控洪水，堤防隱患探測工作已經(jīng)變得十分重要。如何高效、準確地確定隱患的分布位置和大小，已經(jīng)成為堤防隱患探測迫切需要解決的問題 [1]。

20世紀80年代以來，我國開始采用物探方法[2]開展堤防隱患探測的研究工作，如高密度電法[3]、地質(zhì)雷達法[4]、瞬變電磁法[5]等。物探方法是根據(jù)堤防土體與隱患之間的地球物理特性差異來判斷是否存在堤防隱患，在堤防質(zhì)量檢測中起到了重要的作用。

目前，地震CT技術(shù)[6]已經(jīng)被多次應(yīng)用于地球物理勘探的不同工況中[7]。閻生存等[8]對大壩進行ART算法成像研究，推動了CT技術(shù)在大壩隱患檢測中的發(fā)展[9];胡明順等[10]提出井地地震CT探測的方法;李世民等[11]運用地震CT技術(shù)對隧道前方地質(zhì)情況進行超前預(yù)報;靳洪曉等[12]運用地面地震CT技術(shù)對淺層地質(zhì)情況進行勘探研究。但是目前地震CT技術(shù)應(yīng)用在堤防上的案例較少，堤防CT和一般的井間地震CT不同，堤防是一個近似規(guī)則的梯形結(jié)構(gòu)[13]，兩邊均有一定坡度，分為迎水面和背水面。因此，堤防CT模型應(yīng)該是一個梯形模型。

堤防模型的建立，除了考慮近似于梯形結(jié)構(gòu)的地形影響外，還需要考慮堤防的構(gòu)筑材料對波速的影響。堤防主要有土堤、鋼筋混凝土堤、土石混合堤等，其中土堤最為常見，在某些缺乏土料的地區(qū)常采用土石混合堤。堤防的上層一般是不透水覆蓋層，下層一般為透水的砂礫石層和細砂層。堤防的土一般為壤土類土和黏土類土[14]，其中可能夾雜著少量的砂土。這些堤防中的不同構(gòu)筑材料都對波速有著不同程度的影響。

堤防的迎水面和背水面工作條件較好，在采用地震CT進行堤防隱患探測時，容易實現(xiàn)地震波的激發(fā)和接收。本文利用基于MATLAB編寫的程序進行堤防CT的數(shù)值模擬，反演堤防中隱患的位置和大小，從而驗證在堤防上進行地震CT的準確性和可行性。

2 堤防CT基本原理

2.1 離散圖像重建原理

如圖1所示，在Y方向沿梯形左側(cè)（迎水面）安置系列震源、Y方向沿梯形右側(cè)（背水面）安置系列檢波器。

震源產(chǎn)生的地震波被檢波器采集后，以電信號的形式傳輸至主機。通過拾取地震波初至時間，可實現(xiàn)反演這個梯形剖面的速度v（x，y）或慢度f（x，y）=1/v（x，y）。設(shè)第i個地震波的走時為τi，傳播路徑的長度為Li，則

式（1）中df代表弧長微元，速度和傳播路徑都是未知條件，而時間τi是試驗采集到的數(shù)據(jù)，為已知的。在堤防中地下介質(zhì)近似相同，速度場可以看作是不變的，射線的路徑Li看作近似直線。但是，在實際的堤防中地下介質(zhì)情況是不同的，射線路徑也往往是彎曲的。假設(shè)圖像重建區(qū)的網(wǎng)格數(shù)目為M，單元格的慢度依次為f1、f2、…、fM。第i條射線的傳播時間為

通過式（3）計算得出一個速度場的分布矩陣，根據(jù)獲得的數(shù)據(jù)進行反演成像。在計算反演矩陣A時往往是沒有規(guī)則的，且其中的射線只通過了反演區(qū)域的一部分。在實際工程應(yīng)用中常用迭代法求解式（3），重建探測區(qū)域的速度場。

2.2 ART算法

目前，反演的方法很多，主要有以下幾種：矩陣法、付氏變換法、卷積反投影法、代數(shù)重建法（簡稱ART法）[15]。ART算法[16]適用性較廣，直射線和彎曲射線均能使用。

ART算法的基本思想是先給反演區(qū)域中的每個像元都賦予一個初始時間值tij，一般為0;然后將投影結(jié)果沿射線方向反投影回原像元，同時不斷地對反演區(qū)域進行修正;最終在滿足要求時停止修正。假設(shè)fq，ij為第q輪迭代時第i條射線上對于第j像元的慢度（估計值），則第i條射線的理論走時可以表示為

所有射線修正完后，若計算得出走時差的模ΔTq<ε（ε為給定允許誤差），則停止迭代;否則，再進行第q+1輪的迭代，從i=1開始。

2.3 SIRT算法

SIRT算法[17]是由Gilbert首先提出的。該算法與ART算法不同，典型的ART算法用第j條射線的修改值來計算第j+1條射線的修改值，而SIRT算法是把第k輪迭代的所有射線的修改值保存到一個指定空間，在這一輪結(jié)束后求得某種平均修改值Δf-，然后再對每個像元的慢度進行修正，即f（k+1）=f（k）+Δf-，并用于下一輪迭代。

3 軟件介紹

軟件基于MATLAB編寫，反演方法使用ART算法和SIRT算法。MATLAB軟件可以用于矩陣運算、函數(shù)和數(shù)據(jù)繪制、三維成圖和創(chuàng)建圖形用戶界面等[18]，將數(shù)據(jù)進行矩陣化處理，運用MATLAB進行計算，極大地提高了工作效率。

軟件主要分為兩個部分，分別為規(guī)則的梯形堤防模型和有淤背區(qū)的堤防模型。有淤背區(qū)的堤防模型是在規(guī)則梯形軟件的基礎(chǔ)上設(shè)計一個淤背區(qū)。

梯形模型軟件的主要功能是建立一個堤防梯形模型，對此模型進行正演和反演，得出地震波在堤防中的速度分布圖，最后找出堤防隱患的位置，從而提高堤防隱患探測的效率和數(shù)據(jù)解釋的可視化程度。軟件的運行過程如下：第1步，建立一個網(wǎng)格化的梯形模型，在此模型上設(shè)置激發(fā)點和采集點;第2步，分別設(shè)置實際的背景探測速度和假設(shè)的背景探測速度，通過計算得到每條射線經(jīng)過像元的長度，并計算出每條射線的理論時間和實際時間;第3步，將結(jié)果利用ART算法和SIRT算法進行反演，當滿足條件時，停止迭代;第4步，計算結(jié)果成圖，得出探測區(qū)速度分布圖，從圖中的速度分布情況可以看出隱患異常體所在的位置。

4 數(shù)值模擬

堤防觀測系統(tǒng)具備對稱性，因此建立一個梯形堤防數(shù)值模型。異常區(qū)與堤防填筑料的平均速度差對成像質(zhì)量的影響很大，速度差越大，成像質(zhì)量越好。地震波通過堤防時，如果遇裂縫、空洞等異常體，速度會降低。假設(shè)規(guī)則梯形探測區(qū)背景速度為1 000 m/s，模型中異常體像元的速度為900 m/s，如圖2（a）所示，設(shè)計一個規(guī)則的二維梯形模型，網(wǎng)格區(qū)域為5×6，網(wǎng)格劃分為1×1。梯形左側(cè)為迎水面，右側(cè)為背水面。激發(fā)點分布在梯形的左側(cè)，采集點分布在梯形的右側(cè)，設(shè)置激發(fā)點和采集點間距為0.5 m或1 m，重建區(qū)網(wǎng)格數(shù)為26。

淤背區(qū)對堤防工程具有重要的作用，能夠增強堤防的抗?jié)B和抗震穩(wěn)定性。因此，設(shè)計重建區(qū)網(wǎng)格圖形模型為有淤背區(qū)的堤防，右側(cè)有一個高2 m、寬1 m的淤背區(qū)。基于MATLAB使用ART算法和SIRT算法進行反演，最大迭代次數(shù)均為20次。

采用圖像距離測量值來評價圖像重建效果，即歸一化均方根距離測量值d、歸一化平均絕對距離測量值r和最壞情況距離測量值e。d=0表示重建圖像忠實地再現(xiàn)原始圖像，d值愈大表示兩者偏差愈大;r=0說明沒有誤差，r增大，說明誤差增大;e為原始圖像和重建圖像之間的最大平均慢度差。建立3類模型，第1類為含低速體堤防，第2類為層狀堤防，第3類為有淤背區(qū)的層狀堤防，進行對比分析。

4.1 含低速體堤防模型

圖2（a）為激發(fā)點和采集點的間距為1 m的含低速體模型，一共25條射線。圖3（a）為激發(fā)點和采集點的間距為0.5 m的含低速體模型，一共81條射線。

圖2和圖3中的迭代次數(shù)均設(shè)置為20次。激發(fā)點與采集點的間距越小，經(jīng)過每一個像元的射線條數(shù)越多。由表1和圖3可知，在迭代10次后，SIRT算法結(jié)果趨于穩(wěn)定，ART算法的圖像距離測量值還在增大;迭代20次后ART算法的結(jié)果才趨于穩(wěn)定。由此可見，對于堤防含低速體模型，SIRT算法比ART算法收斂速度更快，且收斂性更好。

4.2 層狀堤防模型

層狀堤防ART算法和SIRT反演結(jié)果見圖4。圖4（a）為激發(fā)點和采集點的間距為0.5 m的層狀堤防模型，一共81條射線。探測區(qū)的背景速度為1 000 m/s，模型中像元（層狀）的速度為900 m/s。

由圖4和表2可知，SIRT算法的結(jié)果好于ART算法，且收斂速度更快。ART算法在第10次迭代后，d值還在增大，而SIRT算法d值在逐漸減小。圖4中（c）和（d）迭代次數(shù)均為20次，都可以明顯看出異常體的位置和大小。

4.3 有淤背區(qū)的層狀堤防模型

有淤背的層狀堤防模型反演結(jié)果見圖5和表3。圖5（a）為激發(fā)點和采集點的間距為0.5 m的有淤背區(qū)的層狀堤防模型，共78條射線。探測區(qū)背景速度為1 000 m/s，異常體的速度為900 m/s，在淤背區(qū)有輕微繞射，速度為800 m/s。

由圖5和表3可知，SIRT算法的收斂速度大于ART算法的，且受到淤背區(qū)的干擾較小，上下層都能明顯反演出異常體的位置和大小。

5 結(jié) 語

（1）由于地震波傳播實際路徑并不是一條直線，因此不能完全重建出真實速度場，但根據(jù)速度場相鄰區(qū)域的差異可以識別出異常體在堤防中的空間分布及形態(tài)。

（2）當?shù)螖?shù)相同時，SIRT算法重建圖像的距離測量值比ART算法的小。當距離測量值相近時，SIRT算法重建的迭代次數(shù)少于ART算法的?？梢奡IRT算法比ART算法的迭代收斂性要好，收斂速度快，而且對投影數(shù)據(jù)誤差的敏感度小。

（3）震源和檢波器的排布設(shè)計直接影響射線分布的密度、均勻性，從而影響圖像重建的精度和清晰度。因此，在實際應(yīng)用中，要選擇合適的震源和檢波器的排布，以獲得真實可靠的重建圖像。

利用MATLAB編程軟件進行堤防CT的隱患探測研究，通過數(shù)值模擬，ART算法和SIRT算法反演都得出了較好的結(jié)果，驗證了地震CT進行堤防隱患探測的正確性和可行性。

將地震CT方法應(yīng)用在堤防隱患探測上，提高了隱患探測的效率和精度，并且減少了鉆孔費用。在實際工作中，如果在觀測系統(tǒng)設(shè)計、數(shù)據(jù)采集和數(shù)據(jù)處理三個階段選擇合適的方法，那么堤防CT成像技術(shù)將取得較好的探測結(jié)果。下一步還需要進行實際堤防隱患探測試驗，與實際鉆孔資料對比，進一步驗證堤防CT成像技術(shù)的實用性。

參考文獻：

[1] 周華敏，肖國強，周黎明，等.堤防隱患物探技術(shù)研究現(xiàn)狀與展望[J].長江科學院院報，2019，36（12）：164-168.

[2] 欒明龍.堤防隱患的地球物理模型及其探測與監(jiān)測的思考[J].北京水利，2005（5）：36-38.

[3] 宋朝陽，王銳，李長征，等.高密度電法探測堤防隱患研究[J].人民黃河，2020，42（7）：104-106.

[4] 武桂芝，張寶森，李春江，等.陣列地質(zhì)雷達在黃河堤防隱患探測中的應(yīng)用[J].人民黃河，2020，42（8）：113-116.

[5] 鄧洪亮，謝向文，郭玉松，等.黃河下游堤防工程隱患探測技術(shù)與應(yīng)用[J].地球物理學進展，2008，23（3）：936-941.

[6] 王旭，張晨，由靜，等.地震波CT在隧道地質(zhì)災(zāi)害探測中的應(yīng)用[J].大壩與安全，2010（1）：49-52.

[7] 羅炬，李志海，王海濤.地震層析成像研究方法綜述[J].內(nèi)陸地震，2011，25（4）：311-320.

[8] 閻生存，李珍照，薛桂玉，等.大壩CT中ART算法及其改進探討[J].武漢大學學報（工學版），2001，34（4）：29-34.

[9] 余志雄，薛桂玉，周洪波，等.大壩CT技術(shù)研究概況與進展[J].巖石力學與工程學報，2004，23（8）：1394-1397.

[10] 胡明順，潘冬明，李兵，等.井地地震CT成像數(shù)值模擬研究與應(yīng)用[J].地球物理學進展，2009，24（4）：1399-1404.

[11] 李世民，張曉培，牛建軍，等.地震CT技術(shù)在隧道工程中的應(yīng)用[J].探礦工程（巖土鉆掘工程），2004（9）：63-65.

[12] 靳洪曉，趙永貴，李勤，等.地面地震CT在淺層勘探中的應(yīng)用[J].工程地質(zhì)學報，2000，8（2）：239-243.

[13] 孫東亞，姚秋玲，趙進勇，等.堤防工程建設(shè)技術(shù)進展[J].中國防汛抗旱，2009，19（6）：34-37.

[14] 楊代齊.水利堤防土方填筑方法探討[J].水利技術(shù)監(jiān)督，2015，23（3）：92-95.

[15] 楊文采.地球物理反演的理論與方法[M].北京：地質(zhì)出版社，1997：83-100.

[16] 駱循，朱介壽.ART算法及其在地球物理反演中的應(yīng)用[J].CT理論與應(yīng)用研究，1987（1）：8-16.

[17] 楊利普，徐志萍，李德慶，等.SIRT法層析成像在某核電站巖溶區(qū)的應(yīng)用[J].CT理論與應(yīng)用研究，2014（2）：6-14.

[18] 段成龍，閻長虹，許寶田，等.跨孔地震CT技術(shù)在地鐵工程施工溶洞探測方面的應(yīng)用[J].地質(zhì)論評，2013，59（6）：1242-1248.

【責任編輯 許立新】