鄧 飛 聶煥鑫 沈聯(lián)洪 李 旭 楊 林

(1.成都理工大學計算機與網(wǎng)絡安全學院， 成都 610059； 2.成都交通投資集團有限公司， 成都 610041)

隨著國家經(jīng)濟快速發(fā)展，我國在城市建設方面成效顯著，諸多公路橋梁被不斷建起。但人口數(shù)量的劇增給交通運輸造成了很大壓力，隨著時間推移，在各種因素的影響下，橋梁結構會受到一定損傷。因此橋梁的檢測工作對保障城市的正常交通運輸和推動經(jīng)濟發(fā)展有著重要意義[1-2]。

橋梁檢測作為一門逐漸興起的學科，現(xiàn)有眾多的檢測技術，比如波動層析(CT)檢測技術[3]、光纖傳感檢測技術[4]、探地雷達檢測技術[5]，無人機檢測技術等等[6-7]。其中波動CT檢測技術是利用敲擊產(chǎn)生的聲波或振動波在工程介質(zhì)中的傳播特性，通過檢波器采集穿過橋梁的波并用計算機加以分析，來探查橋梁內(nèi)部結構的一種技術。由于其施工簡便、檢測結果直觀可靠，可以直接反映橋梁內(nèi)部的質(zhì)量、缺陷與損傷等情況，成為了近年來橋梁無損檢測的熱點。

CT技術原為醫(yī)學X射線的理論，后被用于地震層析中檢測地下介質(zhì)的結構情況。同理，也可用于橋梁檢測中。徐蓉等將CT成像技術引入大型橋梁基樁的超聲波檢測中[8]，使檢測結果更直觀、異常圈定更準確；谷昊提出在橋梁檢測中使用聲波CT技術[9]，探討了CT技術對橋梁混凝土檢測的可行性；袁浩等將彈性波層析技術運用于橋梁無損檢測[10]，實踐結果較好地反映出了混凝土澆筑的質(zhì)量情況；Li等研究了一種基于壓縮采樣的超聲波層析技術[11]，進一步減少了測量工作量，可以較快速和精確地檢測橋梁混凝土的缺陷部分。

在橋梁無損檢測中，CT檢測技術有分辨率高、可靠性好等特點，然而以上技術所采用的反演模型大多仍基于矩形網(wǎng)格，不適用于對不規(guī)則橋梁模型的反演計算[12]；且測量代價與工作量等依然較高，在橋梁檢測工程中難以全面實施，無法拾取大量數(shù)據(jù)保證反演精度；此外也未考慮激發(fā)設備無法準確記錄起振時間所造成的時差對小尺度的橋梁反演造成的影響。為此，本文提出一種基于非結構化網(wǎng)格的菲涅爾帶層析反演算法，并結合時差法，通過理論模型與實地測量，驗證了該方法在解決以上提出的橋梁檢測問題上是可行性的。

1 方法原理

橋梁CT檢測是通過激發(fā)振動波穿透混凝土，根據(jù)拾取的振動波走時計算出振動波在橋梁混凝土中的速度分布，從而對混凝土結構進行成像的技術。分為模型正演和層析反演兩個部分，在正演過程中建立橋梁初始模型，使用正演算法計算振動波的理論走時和射線路徑；在反演過程中用實際拾取的振動波初至數(shù)據(jù)與正演結果比較，利用反演算法對理論橋梁模型進行迭代修正，最終得出真實橋梁模型。

基于非結構化網(wǎng)格的時差菲涅爾帶反演方法，是一種使用時差法拾取初至時、在SIRT反演算法上引入菲涅爾帶的的橋梁層析檢測方法。本方法的計算思路如圖1所示，先建立初始橋梁模型和觀測系統(tǒng)，其中橋梁模型的初始速度利用測量數(shù)據(jù)通過式(1)計算得到：

圖1 反演流程Fig.1 Flow chart of tomography

v=(li-lj)/(ti-tj)

(1)

式中：li與lj為單炮中距離起振點最近的兩道；ti與tj為對應的測量初至時。一般認為起振點近道橋梁質(zhì)量良好，因此通過式(1)即可利用不準確的初至時計算出橋梁正常速度，為保證準確性，也可綜合多炮數(shù)據(jù)取均值。

在建立好初始模型后，對模型進行正演計算理論初至；然后利用時差法計算出真實旅行時，將真實旅行時導入橋梁模型，利用基于菲涅爾帶的SIRT反演算法對模型進行反演，并更新模型，重復該過程直至迭代結果滿足終止條件，得到最終反演出的橋梁速度模型。

1.1 基于非結構化網(wǎng)格的正演

在地震層析成像中，常用的正演算法之一為旅行時線性插值法(Linear Travel Interpolation, LTI)，該算法以矩形網(wǎng)格剖分介質(zhì)模型，利用最短路徑算法，通過線性插值公式來向前和向后處理，可以計算出任意網(wǎng)格節(jié)點的旅行時。該方法最初用于地震勘探時，矩形網(wǎng)格對模型的劃分僅在近地表不平處會產(chǎn)生一定的誤差。然而在橋梁檢測中，由于橋梁規(guī)模相對較小、邊緣輪廓不規(guī)則，使用矩形網(wǎng)格劃分會產(chǎn)生較大誤差。

為提高反演結果的整體精確度，本文使用三角網(wǎng)格剖分橋梁模型，對比矩形網(wǎng)格，在相同網(wǎng)格密度下可以更好地貼合橋梁邊緣，因此在正演過程中也需要使用基于非結構化網(wǎng)格的LTI算法[13]，減小了正演射線誤差，配合本文基于三角網(wǎng)格的反演算法，使橋梁層析總體結果的精度得到了提升。反演算法及初至拾取將在1.2、1.3節(jié)中介紹。

1.2 菲涅爾帶反演方法

聯(lián)合迭代重建技術(Simultaneous Iterative Reconstruction Techniques, SIRT)是射線旅行時層析反演的常用方法之一[14]。本文使用的三角網(wǎng)格模型與矩形網(wǎng)格模型在計算上類似，即：

(2)

式中：Sj為第j個三角網(wǎng)格的慢度；q為迭代次數(shù)；Nj為穿過網(wǎng)格j的射線數(shù)目；ΔSjn為穿過網(wǎng)格j的第n條射線對網(wǎng)格j的慢度修改量。

在地震層析反演中，使用了大量檢波設備來采集初至數(shù)據(jù)，大量數(shù)據(jù)的支持使得射線SIRT算法反演出的結果具有較高的精確度，然而在橋梁CT檢測中，受經(jīng)費、施工條件等因素制約，檢波設備不足，導致射線稀疏，模型覆蓋面積小；此外，射線只是地震波的高頻近似，不能反映真實的物理過程，會對反演準確性造成較大影響。

為解決上述問題，本文引入了菲涅爾帶層析反演方法[15]。由于振動波能量在現(xiàn)實中是沿著數(shù)學射線菲涅爾帶范圍內(nèi)的介質(zhì)傳播的，因此在反演過程中加入菲涅爾帶，不僅更能反映振動波的真實傳播情況，也能擴大射線的覆蓋范圍，如圖2所示。

圖2 菲涅爾帶示意Fig.2 The schematic diagram of Fresnel zone

S為炮點，R為檢波點，橢圓條帶為射線SR的菲涅爾帶。A為模型內(nèi)任意點，根據(jù)射線追蹤法可以計算出S到A的旅行時t(S,A)，又根據(jù)振動波傳播的互異性，可計算出A到R的旅行時即為t(R,A)。由此可以得出菲涅爾帶范圍判定式：

|t(S,A)+t(R,A)-t(S,R)|≤1/2f

(3)

式中：t(S,R)即為射線旅行時；f為菲涅爾帶頻率。

根據(jù)式(4)可以計算出模型中任意三角網(wǎng)格是否處于射線菲涅爾帶范圍內(nèi)。

綜上所述，本文提出的三角網(wǎng)菲涅爾帶時差層析反演算法步驟如下：

1)建立橋梁模型，對橋梁模型進行三角剖分，在模型中布設炮點和檢波點；

2)給定初始慢度S0；

3)設定總迭代次數(shù)N，設定反演終止的平均誤差門檻值ξ0；

4)將三角網(wǎng)格的累計慢度修改量設為零；

5)遍歷每個炮點開始反演計算；

6)通過正演計算得到第i個炮點到每個檢波點的旅行時表和射線軌跡，使用時差法對拾取的初至波的旅行時進行計算，得到合成后的旅行時；

7)找出每根射線菲涅爾帶范圍內(nèi)的三角網(wǎng)格，加入表L中；

8)計算表L中每個三角網(wǎng)格的慢度修改量，并累加慢度修改量；

9)重復步驟5)～8)，完成炮點循環(huán)，使用累計慢度修改量修改模型；

11)使用評價函數(shù)：

(4)

12)重復步驟4)～11)，直到滿足ξ<ξ0或迭代次數(shù)大于N，退出迭代，反演計算完成。

對于菲涅爾帶的范圍大小，由于橋梁體積相對較小，使菲涅爾帶頻率f=2 000 Hz進行測量效果較為理想，在上述反演迭代的步驟7)～8)中，需找出射線菲涅爾帶范圍內(nèi)的三角網(wǎng)格，并對其慢度進行修正，具體步驟如下：

1)要找出射線的菲涅爾帶范圍內(nèi)的所有三角網(wǎng)格，先要找出當前射線通過的三角網(wǎng)格，每條射線由多個節(jié)點連接而成，為此先建立一條隊列Q；

2)遍歷射線的每個節(jié)點，得到節(jié)點所在的三角網(wǎng)格，將其加入隊列Q，并標記為菲涅爾帶范圍內(nèi)的三角網(wǎng)格，下述簡稱為加上flag標記，每條射線的flag值不同6，并加入表L當中，這些射線直接通過的三角網(wǎng)格的權值Ω=1；

5)判斷隊列Q是否為空，如果不為空，重復步驟4)，如果隊列Q為空，說明當前射線菲涅爾帶范圍內(nèi)的三角網(wǎng)格已經(jīng)被全部找出并加入表L中，開始下一步計算；

6)遍歷表L中的每個三角網(wǎng)格，將每個網(wǎng)格的權值Ω除以表L中所有三角網(wǎng)格的總權值，得到當前三角網(wǎng)格的權重百分比Ω，則當前三角網(wǎng)格的慢度修改量ΔS=Ω·Δt；

7)表L的所有三角網(wǎng)格計算完畢后，當前射線便計算完成。

1.3 時差法

由于微振層析技術無法準確記錄起振時間[16-17]，單炮記錄的各道初至時間與真實旅行時之間存在相同的毫秒級誤差，導致拾取的初至數(shù)據(jù)不準確。因此本文提出了基于時差的反演方法，在拾取的初至數(shù)據(jù)不準確的情況下，將初至數(shù)據(jù)通過計算化為近似真實的旅行時，可以取得與使用準確數(shù)據(jù)近似的反演效果。

時差法取距離起振點最近的檢波點作為參考，認為該檢波道途經(jīng)的橋梁質(zhì)量良好，其真實旅行時應該與初始模型正演計算出的射線旅行時近似，因此該道在實際中拾取的初至時間與正演旅行時之差即為該炮的時差。單炮第j道的真實旅行時計算如下：

(5)

使用式(5)對實際初至數(shù)據(jù)進行處理后，代入反演計算可以得到較為準確的反演結果。

2 模型試算與實測

2.1 梁柱模型試算

首先在簡單矩形梁柱模型上初步測試菲涅爾帶方法的反演效果，建立10 m×2 m的梁柱初始模型(圖3a)，設置振動波在模型中傳播的速度為3 000 m/s，在梁柱兩邊以1 m為間隔布設檢波點(圖3b中藍點)，并在內(nèi)部設置三列炮點(圖3b中黃點)，先后在高8 m處的梁柱中設置1 m寬、2 000 m/s(圖3c)，1 m寬、1 000 m/s(圖3e)和0.1 m寬、1 000 m/s(圖3g)的低速區(qū)。分別對模型進行正演并將正演結果代入初始模型中，使用菲涅爾帶反演法反演迭代30輪，得到相應反演結果。

根據(jù)圖3d和圖3f可以看出，2 000 m/s和1 000 m/s的1 m低速區(qū)模型反演結果的寬度約為1 m左右，接近于原模型，在區(qū)域邊緣處存在模型1/10左右大小的降速區(qū)。對比圖3h則發(fā)現(xiàn)，當?shù)退賲^(qū)寬度為0.1 m，僅占整體梁柱模型1%的情況下，依然能在8 m高處反演出一片低速區(qū)，寬度在0.1 ～0.3 m之間，但對速度大小的反映略為不準確，此外在區(qū)域邊緣處同樣存在1/10左右的降速區(qū)。據(jù)此可以發(fā)現(xiàn)，隨著低速區(qū)域的縮小，對低速區(qū)的大小和速度的反演精確度逐漸降低，但仍然能夠檢測出速度異常區(qū)域，因此合理設計觀測系統(tǒng)、控制好模型與裂縫的比例，即可對不同大小的裂縫進行精確反演。

a—反演初始模型； b—點位布設；c—1 m寬、2 000 m/s低速區(qū)； d—圖3c的反演結果；e—1 m寬、1 000 m/s低速區(qū)； f—圖3e的反演結果；g—0.1 m寬、1 000 m/s低速區(qū)； h—圖3g的反演結果。圖3 理論梁柱模型試驗結果Fig.3 Test resules of theoretical beam-column model

2.2 橋梁模型試算

下面在橋梁模型上檢驗本文方法的可行性，并與其他方法進行對比。建立最大寬高10 m的橋梁模型，設置振動波在模型中傳播的速度為3 000 m/s，在橋梁邊緣輪廓布設24個檢波點(圖4a中藍點)，并在各處布設炮點(圖4a中黃點)，在橋梁左側和底部設置兩處低速區(qū)，分別為1 000 m/s和500 m/s。對圖4a的模型進行正演，將正演結果代入圖4b的反演初始模型中，分別使用常規(guī)SIRT反演法和菲涅爾帶反演法反演迭代30輪，得到的反演結果如圖4c和圖4d所示，對比兩圖發(fā)現(xiàn)，常規(guī)射線SIRT方法無法正確反演出速度異常區(qū)，而菲涅爾帶SIRT方法的反演結果較準確地揭示出了低速區(qū)的形狀與速度，能夠較好地反映出速度異常區(qū)的形狀和大小。

a—觀測系統(tǒng)模型； b—反演初始模型；c—射線SIRT迭代30輪； d—菲涅爾帶SIRT迭代30輪；e—擾動數(shù)據(jù)迭代30輪； f—時差菲涅爾帶迭代30輪。圖4 理論橋梁模型試驗結果Fig.4 Test results of theoretical bridge model

圖4e為對模型正演旅行時數(shù)據(jù)進行擾動后(分別隨機地對每炮的各道旅行時數(shù)據(jù)進行統(tǒng)一的毫秒級改動)，使用菲涅爾帶方法反演的結果，可以看出反演結果混亂，無法分辨速度異常區(qū)，說明無法準確記錄起振時間對橋梁反演會造成較大影響。對擾動數(shù)據(jù)使用時差法并進行菲涅爾帶反演迭代30輪的結果見圖4f，對比圖4d，發(fā)現(xiàn)使用時差法對擾動數(shù)據(jù)修正后的反演效果與使用準確旅行時數(shù)據(jù)的反演效果近似，說明時差法能夠改善無法準確記錄初至時間的問題。

圖5為反演標準化殘差，由于擾動數(shù)據(jù)的標準化殘差為其余數(shù)據(jù)的十數(shù)倍，且下降速度較慢，因此未添加于圖中?？梢钥吹?，隨著迭代次數(shù)增加，標準化殘差在不斷減小，菲涅爾帶方法的標準化殘差相比SIRT方法減小了1/2左右，而對擾動數(shù)據(jù)使用時差法后，反演標準化殘差也接近于菲涅爾帶法。

圖5 反演標準化殘差Fig.5 Standardized residual of inversion

2.3 實測數(shù)據(jù)處理

理論驗證完畢后，開始實地測量，并對采集到的數(shù)據(jù)進行處理。數(shù)據(jù)采集在一塊帶有裂縫的水泥路面上進行，以裂縫為中心，選取長寬為5 m的區(qū)域。點位布設如圖6所示，在區(qū)域兩側邊緣以0.5 m間隔布設22處檢波點，左下角為1號點，右下角為22號點。區(qū)域內(nèi)以1 m為間隔設置4列共24處炮點(區(qū)域內(nèi)標記處)，在炮點處用鐵錘砸擊鐵餅，各檢波點即能接受到振動波并將數(shù)據(jù)反饋至工程探測儀。

圖6 現(xiàn)場施工照片F(xiàn)ig.6 The photo of site operation

圖7為采集的振動波初至波形，從近炮道開始，拾取近炮道波形頂部對應的時間，并依次拾取后續(xù)道的初至(圖7中紅線)，即得到當前炮各道對應的初至。

a—第1炮記錄； b—第24炮記錄。圖7 初至波波形Fig.7 Oscillogram of first arrival

初至拾取完畢后，建立如圖8a 所示的初始模型。利用式(1)對水泥路面的速度計算得出大致速度為2 650 m/s，并寫入初始模型中。隨后使用時差法對初至數(shù)據(jù)進行處理，并代入初始模型中進行反演計算。

a—觀測系統(tǒng)反演初始模型； b—射線SIRT反演效果；c—菲涅爾帶SIRT反演效果； d—時差菲涅爾帶反演效果。圖8 實測反演結果Fig.8 Inversion results of actual measurement

圖8b為使用拾取的原始初至數(shù)據(jù)，并使用常規(guī)SIRT算法對初始模型反演迭代20次的結果，圖8c和圖8d分別為菲涅爾帶方法和時差菲涅爾帶方法迭代20次的結果?？梢钥闯觯河捎谏渚€不夠密集，圖8b中SIRT方法對模型速度的修改主要集中于炮點及檢波點附近，無法辨認速度異常區(qū)的形狀。根據(jù)圖8c可以看出，在模型中央存在帶狀低速區(qū)，但由于拾取的數(shù)據(jù)存在誤差，模型反演結果的左上部存在明顯的大范圍高速區(qū)及小部分低速區(qū)，與測量地實際情況不符。而反觀圖8d，在對拾取數(shù)據(jù)使用時差法處理后，使用菲涅爾帶方法反演的結果未對其余區(qū)域過度修正，較準確地反映出了模型中央的裂縫。

圖9為實測數(shù)據(jù)反演迭代過程中的標準化殘差，可以看出：在不使用時差法的情況下，反演標準化殘差相對較大，菲涅爾帶方法的殘差曲線優(yōu)于常規(guī)SIRT方法，而時差菲涅爾帶法的標準化殘差則遠小于以上二者，并且從圖8的反演效果也可以看出，時差菲涅爾帶法在對實際模型的修正上同樣能取得良好效果。

圖9 實測標準化殘差Fig.9 Standardized residual of actual measurement

3 結 論

利用基于三角網(wǎng)的時差菲涅爾帶反演方法，分別對理論數(shù)據(jù)和實際測量數(shù)據(jù)進行了反演計算，通過對比常規(guī)SIRT方法、菲涅爾帶方法以及時差菲涅爾帶方法的反演結果，得到以下結論：

1)針對橋梁檢測提出的基于三角網(wǎng)的時差菲涅爾帶反演方法，通過對模擬梁柱模型、橋梁模型和實際測量數(shù)據(jù)的計算，證明了該方法在建筑結構梁、柱與橋梁檢測領域的可行性。

2)理論試驗和實際測量結果表明，菲涅爾帶方法在檢波點較少的情況下也能夠?qū)δＰ瓦M行較正確的反演，而時差法則可以有效減小初至時間記錄不準確造成的誤差。基于三角網(wǎng)的時差菲涅爾帶反演方法，對比常規(guī)射線SIRT方法，在對模型形狀和速度的反演效果上有明顯改善，可以應用于實際橋梁檢測工程中。