橋梁檢測(cè)中非結(jié)構(gòu)化網(wǎng)格層析正演算法研究

李 旭， 房 鵬， 祁 勛， 楊 林， 趙 宇， 付 濤， 黃玲菱

(1.成都交通投資集團(tuán)有限公司，成都 610041；2.成都理工大學(xué) 信息科學(xué)與技術(shù)學(xué)院，成都 610059)

0 引言

隨著國(guó)民經(jīng)濟(jì)的飛速發(fā)展，我國(guó)的交通運(yùn)輸基礎(chǔ)設(shè)施建設(shè)取得了巨大的成就。特別是橋梁工程的試驗(yàn)檢測(cè)工作是質(zhì)量監(jiān)督和管理保障中核心所在。橋梁檢測(cè)技術(shù)是一門正在發(fā)展的新興學(xué)科，近幾年各類技術(shù)手段層出不窮，如機(jī)械檢測(cè)技術(shù)[1]、振動(dòng)檢測(cè)技術(shù)[2]、電磁檢測(cè)技術(shù)[3]、超聲波檢測(cè)技術(shù)[4]、攝像檢測(cè)技術(shù)[5]等。

振動(dòng)檢測(cè)技術(shù)是指研究機(jī)械振動(dòng)或由機(jī)械振動(dòng)引起的應(yīng)力波在工程介質(zhì)中傳播特性，通過(guò)科學(xué)分析獲取橋梁的技術(shù)參數(shù)和計(jì)量信息的一種新興技術(shù)[6-8]。基于振動(dòng)和沖擊原理的振動(dòng)檢測(cè)技術(shù)，已成為近年來(lái)橋梁動(dòng)態(tài)無(wú)損檢測(cè)的一個(gè)熱點(diǎn)。其中混凝土橋梁振動(dòng)波層析成像檢測(cè)技術(shù)由于其施工快捷、分辨率高、可靠性好、圖像直觀以及可直接反應(yīng)混凝土力學(xué)性狀等等優(yōu)點(diǎn)，尤其得到檢測(cè)機(jī)構(gòu)的重視。該技術(shù)通過(guò)波的傳播特性來(lái)評(píng)價(jià)混凝土的強(qiáng)度、密實(shí)性、均勻性、內(nèi)部缺陷與損傷等，可用于橋梁施工中的質(zhì)量控制、運(yùn)營(yíng)橋梁健康與病害檢測(cè)以及加固處理效果的檢測(cè)評(píng)價(jià)等。

波動(dòng)層析成像檢測(cè)技術(shù)中波動(dòng)層析正演是其核心技術(shù)之一。Moser[9]對(duì)射線追蹤法進(jìn)行全面的研究；國(guó)內(nèi)劉洪等[10]提出了計(jì)算最小走時(shí)的界面網(wǎng)格方法；張建中[11]等研究了射線最短路徑追蹤算法的正確性與效率；魯彬[12]等利用迭代法優(yōu)化了網(wǎng)絡(luò)最短路徑射線算法，計(jì)算效率得到很大程度的提高；劉潤(rùn)澤等[13]提出了小尺度下的三角網(wǎng)聲波層析成像技術(shù)。以上技術(shù)所采用的參數(shù)化模型大多基于矩形網(wǎng)或計(jì)算效率較低，但是由于橋梁整體為不規(guī)則模型，常規(guī)的矩形網(wǎng)格技術(shù)已不適用于波動(dòng)的正演計(jì)算，這里提出一種自適應(yīng)非結(jié)構(gòu)化網(wǎng)格波動(dòng)層析正演算法，并結(jié)合理論模型驗(yàn)證算法的正確性與計(jì)算精度。

1 方法原理

基于線性插值的最小旅行時(shí)正演算法，是地震勘探層析正演中的常用算法，該算法將地下介質(zhì)劃分成離散的矩形網(wǎng)格單元，利用最短路徑算法從震源點(diǎn)開(kāi)始逐層向外尋求子震源點(diǎn)，并且利用LTI(Linear Travel Time Interpolation Ray Tracing)插值計(jì)算擴(kuò)展網(wǎng)格結(jié)點(diǎn)的旅行時(shí)。在地震勘探中，由于地下介質(zhì)是無(wú)限半空間，因此利用矩形網(wǎng)格僅在起伏劇烈的地表會(huì)產(chǎn)生一定的誤差。然而在橋梁檢測(cè)中，由于橋梁模型外形輪廓復(fù)雜如果使用矩形網(wǎng)格來(lái)近似將會(huì)產(chǎn)生較大誤差，為此筆者提出一種自適應(yīng)非結(jié)構(gòu)化網(wǎng)格的層析正演算法，算法基于非結(jié)構(gòu)化的三角網(wǎng)格模型，在不改變格網(wǎng)密度的前提下，三角網(wǎng)格可以表現(xiàn)復(fù)雜的模型結(jié)構(gòu)，通過(guò)測(cè)試對(duì)比，三角網(wǎng)相比矩形網(wǎng)格，既減少了規(guī)則格網(wǎng)帶來(lái)的數(shù)據(jù)冗余，同時(shí)在計(jì)算效率方面又優(yōu)于矩形網(wǎng)格模型。

與矩形網(wǎng)格LTI算法類似，基于非結(jié)構(gòu)化的三角網(wǎng)格LTI算法的核心也是使用旅行時(shí)線性插值方法，用已知最小旅行時(shí)的兩個(gè)頂點(diǎn)估算待插值點(diǎn)處的最小旅行時(shí)，其原理如圖1所示。已知A點(diǎn)和B點(diǎn)的坐標(biāo)(xA,zA)、(xB,zB)和各自旅行時(shí)tA、tB，則線段AB之間任意一點(diǎn)D(xD,zD)的旅行時(shí)tD可由A、B兩點(diǎn)的旅行時(shí)tA、tB線性插值得到，即

(1)

若射線從A，B之間的D點(diǎn)穿過(guò)到達(dá)C點(diǎn)，則C點(diǎn)的旅行時(shí)tC可表示為：

(2)

其中：v為介質(zhì)的波傳播速度；φ為線段AB與AC之間的夾角。

式(2)中C點(diǎn)的旅行時(shí)是關(guān)于r的函數(shù)。根據(jù)費(fèi)馬原理，射線從D點(diǎn)傳到C點(diǎn)應(yīng)滿足:

(3)

聯(lián)立公式(2)、式(3)可以解出:

(4)

圖1 非規(guī)則三角網(wǎng)自動(dòng)剖分

將式(4)代入式(2)可得C點(diǎn)最小旅行時(shí)為:

(5)

基于非結(jié)構(gòu)化的三角網(wǎng)格LTI算法基本的思路為：①?gòu)恼鹪袋c(diǎn)出發(fā)，計(jì)算出包含震源點(diǎn)三角形頂點(diǎn)的旅行時(shí)；②以三角形的每條邊通過(guò)廣度優(yōu)先搜索逐層向外擴(kuò)展插值計(jì)算，尋找旅行時(shí)最小頂點(diǎn)；③以這些頂點(diǎn)作為子震源點(diǎn)，重復(fù)之前的搜索策略，最終得到所有三角網(wǎng)格頂點(diǎn)上的旅行時(shí)。為了便于介紹詳細(xì)的算法過(guò)程，這里先約定如下專用名詞：

1)可到達(dá)點(diǎn)：已經(jīng)計(jì)算過(guò)旅行時(shí)的網(wǎng)格結(jié)點(diǎn)，但不確定求得了最小旅行時(shí)。

2)固定點(diǎn)：已經(jīng)確定求得了最小旅行時(shí)的網(wǎng)格結(jié)點(diǎn)。

3)子震源點(diǎn)：在每一輪計(jì)算中，從當(dāng)前所有可到達(dá)的結(jié)點(diǎn)中找出的，未做過(guò)子震源且旅行時(shí)最小的結(jié)點(diǎn)。

4)擴(kuò)展邊表：由子震源和另一個(gè)固定點(diǎn)構(gòu)成的邊，可用于插值計(jì)算該邊對(duì)點(diǎn)的旅行時(shí)。

根據(jù)上述名詞解釋，提出新一代自適應(yīng)三角網(wǎng)格LTI的算法，其計(jì)算流程如下：

1)計(jì)算包含震源點(diǎn)的三角形三個(gè)頂點(diǎn)的旅行。將三個(gè)頂點(diǎn)標(biāo)記為固定點(diǎn)，并且將它們加入到可到達(dá)點(diǎn)表Q。

2)判斷表Q是否為空。若為空，算法結(jié)束；否則進(jìn)行步驟3)。

3)從Q中找出旅行時(shí)最小的點(diǎn)Pi，若Pi不是固定點(diǎn)，將其標(biāo)記為固定點(diǎn)。

4)以Pi為子震源點(diǎn)，遍歷所有與Pi點(diǎn)共邊的鄰接點(diǎn)Pj。

圖2 LTI算法流程圖

圖3 邊長(zhǎng)0.2 m網(wǎng)格射線追蹤結(jié)果

5)若Pj是固定點(diǎn)，且ΔPiPjPk的另一點(diǎn)Pk不是固定點(diǎn)，則利用線性旅行時(shí)插值公式(LTI)，以Pi—Pj為擴(kuò)展邊，插值計(jì)算Pk點(diǎn)的旅行時(shí)tk。

6)如果tk小于Pk點(diǎn)已有旅行時(shí)Tk，則令Tk=tk。

7)如果Pk原來(lái)不是可達(dá)點(diǎn)，則將Pk加入到Q表中。

8)將Pi點(diǎn)從表Q中去除，返回步驟2)。

綜上步驟的流程見(jiàn)圖2。

2 模型試算

試算1采用相同射線追蹤算法分別對(duì)自適應(yīng)三角網(wǎng)和規(guī)則矩形網(wǎng)進(jìn)行正演模擬，對(duì)比兩種網(wǎng)格的正演結(jié)果，驗(yàn)證算法的正確性與計(jì)算精度。

模型1采用邊長(zhǎng)0.2 m網(wǎng)格，正中心震源激發(fā)，其射線路徑與分布如圖3所示，兩種網(wǎng)格的射線旅行時(shí)誤差分析如圖4所示。圖3(a)為三角網(wǎng)的射線分布，圖3(b)為矩形網(wǎng)的射線分布，兩種網(wǎng)格分布基本一致，但是三角網(wǎng)的射線更加均勻，矩形網(wǎng)在右上角接近45°處出現(xiàn)異常區(qū)域。圖4中可射線旅行時(shí)誤差曲線可反映出這種情況。三角網(wǎng)由于射線分布均勻，累積誤差為1.057 36E-06，其誤差曲線趨于平緩且穩(wěn)定，便于數(shù)值分析與校正。矩形網(wǎng)格累積誤差為1.601 45E-06，略差于三角網(wǎng)的累積誤差，其誤差曲線起伏較大，不利于下一步的旅行時(shí)校正。

圖4 邊長(zhǎng)0.2 m網(wǎng)格射線旅行時(shí)誤差

模型2采用邊長(zhǎng)0.1 m網(wǎng)格，正中心震源激發(fā)，其射線路徑與分布如圖5所示，兩種網(wǎng)格的射線旅行時(shí)誤差分析如圖6所示。圖5(a)為三角網(wǎng)的射線分布，圖5(b)為矩形網(wǎng)的射線分布，與0.2 m網(wǎng)格分布基本一致。圖6中三角網(wǎng)射線旅行時(shí)累積誤差為5.455 15E-07。矩形網(wǎng)格射線旅行時(shí)累積誤差為6.027 42E-07。對(duì)比0.2 m網(wǎng)格結(jié)果可看出兩種網(wǎng)格其累計(jì)誤差均降低了一個(gè)量級(jí)，但是三角網(wǎng)射線旅行時(shí)誤差依舊小于矩形網(wǎng)。

圖5 邊長(zhǎng)0.1 m網(wǎng)格射線追蹤結(jié)果

圖6 邊長(zhǎng)0.1 m網(wǎng)格射線旅行時(shí)誤差

試算2模擬不規(guī)則橋墩上的彈性波正演，分別進(jìn)行均速和變速模型的正演試算。圖7、圖8分別為兩種橋墩模型的射線路徑。其中圖7中彈性波波速為3 000 m/s，藍(lán)色點(diǎn)為傳感器位置，黃色點(diǎn)為震源位于(5, 5)處，射線路徑?jīng)]有出現(xiàn)由于波速變化導(dǎo)致的彎曲現(xiàn)象，由于射線折射點(diǎn)存在，邊界上的傳感器都能接收到震動(dòng)信號(hào)。圖8為變速橋墩模型，其速度變化較大，圖8中藍(lán)色點(diǎn)為傳感器位置，黃色點(diǎn)為震源位于(5, 5)處，射線路徑出現(xiàn)由于波速變化導(dǎo)致的射線彎曲現(xiàn)象，同時(shí)由于折射點(diǎn)存在，邊界上的傳感器皆能收到震動(dòng)信號(hào)。

圖7 勻速橋墩模型射線追蹤效果

圖8 變速橋墩模型射線追蹤效果

3 結(jié)論

筆者利用自適應(yīng)非結(jié)構(gòu)化網(wǎng)格波動(dòng)層析正演模擬橋墩內(nèi)部彈性波的傳播，建立一種非結(jié)構(gòu)化三角網(wǎng)的射線正演模型，通過(guò)規(guī)則模型和非規(guī)則模型試算，分別與矩形網(wǎng)格進(jìn)行正確性和計(jì)算精度對(duì)比,得出以下結(jié)論：

1)當(dāng)采用規(guī)則均勻速度模型檢測(cè)方式時(shí)，在正演劃分的水平和豎直網(wǎng)格數(shù)相等的情況下，自適應(yīng)三角網(wǎng)分割方法從射線分布特征、誤差曲線走勢(shì)分析、累積誤差來(lái)看是優(yōu)于矩形網(wǎng)格。兩種方法的網(wǎng)格劃分越小，計(jì)算精度越高，累積誤差量越低。三角網(wǎng)的誤差曲線趨勢(shì)較矩形網(wǎng)格較為平緩，累積誤差量較小，利于進(jìn)一步數(shù)值分析與校正。

2)當(dāng)采用非規(guī)則均速與變速度模型檢測(cè)方式時(shí)，在正演劃分的水平和豎直網(wǎng)格數(shù)相等的情況下，本文提出的三角網(wǎng)格LTI算法都能夠有效的計(jì)算射線路徑，同時(shí)反映射線在速度變化中模型出現(xiàn)的彎曲現(xiàn)象。表明通過(guò)自適應(yīng)三角網(wǎng)的劃分，結(jié)合LTI射線追蹤算法對(duì)橋梁內(nèi)部缺陷的檢測(cè)是可靠的。