張 彧，程 華，黃曉寒，袁 野

(后勤工程學(xué)院 軍事土木工程系， 重慶 401311)

基于陣列傳感器相關(guān)性的CT圖像重建研究

張 彧，程 華，黃曉寒，袁 野

(后勤工程學(xué)院 軍事土木工程系， 重慶 401311)

彈性波CT層析成像技術(shù)可以在不損傷混凝土構(gòu)件內(nèi)部結(jié)構(gòu)的情況下對(duì)其內(nèi)部缺陷的信息進(jìn)行提取，再利用相關(guān)算法重建圖像。對(duì)接收點(diǎn)陣列間的相關(guān)性加以考慮，通過(guò)線(xiàn)性變換的方法組合方程，進(jìn)而得到了一種重構(gòu)方程的反演算法。由此將線(xiàn)性靜定方程轉(zhuǎn)變?yōu)槌o定方程，減小系數(shù)矩陣的條件數(shù)，從而達(dá)到了減小測(cè)量誤差抑制噪聲影響的目的。數(shù)值模擬結(jié)果表明:在相同條件下，反演得到的圖像更為清晰。從實(shí)驗(yàn)數(shù)據(jù)上看，重建所得的圖像噪聲更小，對(duì)缺陷判別的準(zhǔn)確度更高。

彈性波；反演算法；方程重構(gòu)；圖像重建

CT層析成像技術(shù)是指通過(guò)物體外部檢測(cè)到的數(shù)據(jù)重建物體內(nèi)部(橫截面)信息的技術(shù)，它把不可分割的對(duì)象假想地切成一系列薄片，分別給出每一片上的物體的圖像，然后把這一系列圖像疊加起來(lái)，就得到物體內(nèi)部的圖像。在CT層析成像技術(shù)所建立的走時(shí)方程中，系數(shù)矩陣為大型稀疏矩陣，具有較高的病態(tài)性，方程組的解受誤差的影響較大，所以可以通過(guò)一維傳感器陣列測(cè)點(diǎn)之間的相關(guān)性關(guān)系重構(gòu)彈性波CT方程，將線(xiàn)性靜定方程組轉(zhuǎn)變?yōu)槌o定方程組求解，從而使方程組的解具有較好的精度和穩(wěn)定性。

在彈性波CT方程中只利用了激勵(lì)點(diǎn)與傳感器測(cè)點(diǎn)的相關(guān)性，忽視了傳感器陣列測(cè)點(diǎn)間的相關(guān)性。傳感器陣列測(cè)點(diǎn)的相關(guān)性可通過(guò)本文提出的基于多時(shí)窗能量比和廣義二次互相關(guān)的時(shí)延估計(jì)得到。首先識(shí)別信號(hào)起跳點(diǎn)的大致范圍，利用多時(shí)窗能量比方法對(duì)信號(hào)的起跳點(diǎn)進(jìn)行估計(jì)，然后將預(yù)估的起跳點(diǎn)前后一定范圍以外的信號(hào)的幅值全部置零，再利用廣義二次互相關(guān)對(duì)處理后的信號(hào)進(jìn)行時(shí)延估計(jì)，最后再結(jié)合希爾伯特變換對(duì)廣義二次互相關(guān)的峰值進(jìn)行銳化處理得到時(shí)延。該方法能快速準(zhǔn)確地對(duì)傳感器陣列測(cè)點(diǎn)相關(guān)性進(jìn)行估計(jì)，具有較高的穩(wěn)定性，同時(shí)在較低信噪比的情況下也能獲得比較理想的結(jié)果，為彈性波CT方程反演成像提供了良好的基礎(chǔ)。

一維傳感器陣列相關(guān)性示意圖如圖1所示。該方法通過(guò)利用測(cè)點(diǎn)間的相關(guān)性即測(cè)點(diǎn)間的時(shí)延來(lái)達(dá)到降低噪聲和誤差帶來(lái)的影響，其中加速度傳感器采集得到的信號(hào)1和信號(hào)2如圖2所示。

圖1 一維傳感器陣列相關(guān)性示意圖

彈性波時(shí)延提取技術(shù)對(duì)于研究所測(cè)物體內(nèi)部特性來(lái)說(shuō)是一項(xiàng)非常重要的技術(shù)，且由于初至波形變化較大，相鄰道的波形互相干擾，給時(shí)延的判斷帶來(lái)了巨大困難。目前，已有許多學(xué)者針對(duì)此問(wèn)題提出了相關(guān)的解決辦法。例如徐鈺等[1]在能量比法的基礎(chǔ)上提出了多時(shí)窗能量比初至拾取算法，該算法具有較高的穩(wěn)定性。左國(guó)平等[2]利用時(shí)窗滾動(dòng)來(lái)計(jì)算能量比值，通過(guò)尋找最大的能量比值實(shí)現(xiàn)對(duì)時(shí)延的估計(jì)。向陽(yáng)等[3]將功率譜的相位譜應(yīng)用于時(shí)間延遲估計(jì)。張軍華等[4]將小波變換和能量比方法結(jié)合進(jìn)行初至走時(shí)的提取。王春艷等[5]針對(duì)多通道陣列聲波信號(hào)時(shí)延估計(jì)的問(wèn)題提出了一種平均廣義互相關(guān)函數(shù)算法。針對(duì)這些算法中存在穩(wěn)定性不好、精度不高等問(wèn)題，本文通過(guò)結(jié)合已有的時(shí)延估計(jì)方法和處理手段，成功實(shí)現(xiàn)了對(duì)時(shí)延比較準(zhǔn)確的估計(jì)。

圖2 信號(hào)相關(guān)性示意圖

1 多時(shí)窗能量比法

能量比值法由Coppens[6]提出，其定義為一周期內(nèi)的信號(hào)能量與總時(shí)窗能量的比值，即

(1)

式中：R(τ)為能量比值函數(shù)；x(t)為實(shí)際信號(hào)記錄的振幅；L為視周期的長(zhǎng)度。

能量比值法對(duì)在初至波形變化較小區(qū)域時(shí)可以取得較好的效果，但該方法對(duì)于特征變化明顯的初至波形的時(shí)延估計(jì)效果較差。針對(duì)此方法存在的不足，有學(xué)者提出采用滑動(dòng)時(shí)窗能量比法，即定義如圖3所示的前后時(shí)窗的能量比值。

圖3 滑動(dòng)時(shí)窗能量比法示意圖

在實(shí)際工程應(yīng)用中，信號(hào)的采集具有固定采樣頻率，所以對(duì)式(1)做離散化處理后得到

(2)

式中：T1為第一時(shí)窗起點(diǎn)；T0為第1時(shí)窗終點(diǎn)、第2時(shí)窗起點(diǎn)；T2為第2時(shí)窗終點(diǎn)。由于起跳點(diǎn)附近存在少數(shù)樣點(diǎn)的幅值接近于0，使得式(2)分母趨于0，為避免此情形的發(fā)生，對(duì)式(2)分子分母同時(shí)加上穩(wěn)定因子，得

(3)

式中：A為該信號(hào)道的相對(duì)能量；ω為該信號(hào)道的樣點(diǎn)總數(shù)；α為穩(wěn)定系數(shù)，通常取0.5～2.0。

滑動(dòng)時(shí)窗能量比法僅僅是通過(guò)尋找能量比的最大值以獲得初至走時(shí)，這種方法存在固有缺陷，因?yàn)槌踔敛ㄌ幍哪芰勘炔灰欢ㄟ_(dá)到最大，從而導(dǎo)致初至走時(shí)的拾取錯(cuò)誤。因此，本文采用多時(shí)窗能量比法時(shí)不僅考慮了能量比極大值，同時(shí)也考慮了能量比次級(jí)值，對(duì)于低信噪比信號(hào)初至走時(shí)的拾取精度較好，穩(wěn)定性也較好。

多時(shí)窗能量比法采用4個(gè)滑動(dòng)時(shí)窗，如圖4所示，其中L1、L2、L3、L4分別為第1～4時(shí)窗。

第1時(shí)窗與第2時(shí)窗能量比為

(4)

第1時(shí)窗與第3時(shí)窗能量比為

(5)

第4時(shí)窗與第2時(shí)窗能量比為

(6)

式中：t為當(dāng)前采樣點(diǎn)；d為后時(shí)窗與第3時(shí)窗的間隔；α為穩(wěn)定系數(shù)，一般取值為0.5～2.0。

圖4 多時(shí)窗能量比法示意圖

計(jì)算步驟為：

1) 根據(jù)式(4)計(jì)算A1(t)，搜索A1(t)極大值A(chǔ)max，并記錄其對(duì)應(yīng)的時(shí)間T1。

圖5為彈性波響應(yīng)信號(hào)，圖6為根據(jù)多時(shí)窗能量比法拾取的起跳點(diǎn)。對(duì)求得初至走時(shí)的信號(hào)進(jìn)行處理，保留初至走時(shí)點(diǎn)附近一定范圍內(nèi)的信號(hào)，范圍外的信號(hào)幅值全部置零后得到響應(yīng)信號(hào)，結(jié)果見(jiàn)圖7。通過(guò)以上的處理能夠大致確定起跳點(diǎn)的范圍，同時(shí)對(duì)冗余的信號(hào)成分進(jìn)行置零替換，可大幅度提升對(duì)于初至走時(shí)的判斷精度。

圖5 彈性波響應(yīng)信號(hào)

圖6 多時(shí)窗能量比法拾取起跳點(diǎn)

圖7 處理后信號(hào)

2 廣義二次互相關(guān)的希爾伯特變換

2.1 二次互相關(guān)

假設(shè)信號(hào)模型為：

(7)

其中：s(n)為彈性波信號(hào)；D為延遲時(shí)間；n1(n)、n2(n)為加性噪聲。

x1(n)和x2(n)的相關(guān)函數(shù)表示為

R12(τ)=E[x1(n)x2(n-τ)]=

E[s(n)s(n-τ-D)]+E[s(n)n2(n-τ)]+

E[s(n-τ-D)n1(n)]+

E[n1(n)n2(n-τ)]

(8)

假設(shè)信號(hào)與噪聲以及噪聲與噪聲之間都不相關(guān)，則式(8)變?yōu)?/p>

R12(τ)=E[s(n)s(n-τ-D)]=Rss(τ-D)

(9)

對(duì)于式(8)，當(dāng)τ-D=0時(shí)，R12(τ)取得最大值，此時(shí)相關(guān)函數(shù)的峰值點(diǎn)對(duì)應(yīng)的橫坐標(biāo)就是時(shí)間延遲。當(dāng)信號(hào)信噪比較低、與噪聲相關(guān)時(shí)，利用一次互相關(guān)無(wú)法對(duì)時(shí)延進(jìn)行估計(jì)。針對(duì)此問(wèn)題，考慮使用二次相關(guān)法，將x1(n)的自相關(guān)函數(shù)R11(τ)和x1(n)及x2(n)的互相關(guān)函數(shù)R12(τ)再做相關(guān)，從而提高信號(hào)信噪比與分辨力，R11和R12仍然是時(shí)間的函數(shù)，則其相關(guān)函數(shù)為

RRR(τ)=E[R11(n)R12(n-τ)]=

E[RSS(n)+Rsn1(n)+Rn1s(n)+Rn1n1(n)]·

[RSS(n-D+τ)+Rn1s(n-D+τ)+

Rn1s(n-D+τ)+Rsn2(n+τ)+

Rn1n2(n+τ)]

(10)

彈性波信號(hào)與噪聲的相關(guān)函數(shù)可以忽略，則式(8)變?yōu)?/p>

RRR(τ)=RRS(τ-D)+RRN(τ)

(11)

式中：RRS為純信號(hào)做二次相關(guān)；RRN為噪聲做二次自相關(guān)。對(duì)于式(9)，當(dāng)τ=D時(shí)，RRR(τ)取得最大值，其峰值點(diǎn)對(duì)應(yīng)的橫坐標(biāo)為時(shí)間延遲。二次相關(guān)法的優(yōu)勢(shì)在于能夠進(jìn)一步降低噪聲對(duì)信號(hào)的影響，可以在更低信噪比環(huán)境下估計(jì)時(shí)間延遲。

2.2 廣義互相關(guān)

廣義互相關(guān)算法[7-8]是利用兩個(gè)信號(hào)的廣義相關(guān)函數(shù)來(lái)估計(jì)時(shí)延。為了獲得更好的時(shí)延估計(jì)精度，先利用頻域加權(quán)函數(shù)對(duì)信號(hào)進(jìn)行預(yù)濾波，然后再進(jìn)行相關(guān)性的計(jì)算，檢測(cè)出相關(guān)函數(shù)峰值的橫坐標(biāo)即為時(shí)延。定義經(jīng)過(guò)加權(quán)的廣義相關(guān)函數(shù)為

(12)

式中：Gx1x2(ω)=E[x1(ω)×x2(ω)]為信號(hào)x1和x2的互功率譜；ψ12(ω)為廣義加權(quán)函數(shù)。本文采用的廣義加權(quán)函數(shù)為平滑相干變換窗：

(13)

其中，Gx1(ω)、Gx2(ω)分別為信號(hào)x1和x2的功率譜。那么廣義二次相關(guān)函數(shù)的算法流程如圖8所示。

圖8 廣義二次互相關(guān)算法流程

2.3 希爾伯特變換

希爾波特變換定義如下

(14)

希爾波特變換[9-10]是把相關(guān)函數(shù)的偶對(duì)稱(chēng)性轉(zhuǎn)換成奇對(duì)稱(chēng)性，將相關(guān)法的峰值檢測(cè)轉(zhuǎn)換成過(guò)零檢測(cè)。希爾波特變換中的過(guò)零點(diǎn)的位置正好對(duì)應(yīng)相關(guān)算法中的相關(guān)峰值的位置。通過(guò)將相關(guān)函數(shù)與其希爾伯特變換的絕對(duì)值進(jìn)行相減，既保留了相關(guān)峰值點(diǎn)，又使得峰值附近的相關(guān)值減小，達(dá)到了銳化峰值點(diǎn)的作用，確保了時(shí)延估計(jì)的精度。廣義二次互相關(guān)的希爾伯特變換示意圖如圖9所示。

圖9 希爾伯特變換示意圖

2.4 方程重構(gòu)的實(shí)現(xiàn)與求解

首先對(duì)于選取哪種形式的線(xiàn)性變換方程作為重構(gòu)的方程進(jìn)行了討論，考慮到在獲取走時(shí)兩信號(hào)時(shí)延越大越容易獲得精確的估計(jì)，而且在產(chǎn)生相同的時(shí)延估計(jì)誤差的情況下，對(duì)于時(shí)延越大的兩信號(hào)其誤差所占的比例越小，對(duì)整體的結(jié)果影響較小。所以在同一激勵(lì)條件下，本文選取兩測(cè)點(diǎn)間時(shí)延較大的進(jìn)行線(xiàn)性變換作為附加的方程。如圖10所示，選取接收點(diǎn)1、8的走時(shí)差作為相差，以此類(lèi)推，補(bǔ)充方程的形式如下：

ΔDRm-RnS=ΔTmaxRm-Rn

(15)

式中：ΔDRm-Rn表示相同激勵(lì)下第m個(gè)接收點(diǎn)和第n個(gè)接收點(diǎn)的系數(shù)差；S表示慢度向量； ΔTmaxRm-Rn表示相同激勵(lì)下第m個(gè)接收點(diǎn)和第n個(gè)接收點(diǎn)相關(guān)函數(shù)的極大值。

對(duì)添加方程的數(shù)量進(jìn)行了討論：實(shí)測(cè)過(guò)程當(dāng)中，方程數(shù)量的增加會(huì)大大增加計(jì)算量并帶來(lái)誤差，進(jìn)而影響結(jié)果的精度，所以對(duì)于8×8的網(wǎng)格單元，附加8個(gè)方程進(jìn)行分析，補(bǔ)充后的方程表達(dá)式如下：

圖10 一維傳感器陣列相關(guān)性示意圖

3 數(shù)值算例

3.1 算例1

如圖11所示，設(shè)計(jì)2個(gè)模型截面尺寸為800 mm×800 mm，缺陷尺寸為200 mm×200 mm、100 mm×100 mm的單缺陷數(shù)值模型，并將其劃分為8×8的方形網(wǎng)格。通過(guò)計(jì)算均值、標(biāo)準(zhǔn)差、相對(duì)誤差進(jìn)行對(duì)比，其對(duì)比結(jié)果分別如表1、2所示。

圖11 不同大小缺陷數(shù)值模型

經(jīng)對(duì)比發(fā)現(xiàn)：對(duì)于不同尺寸的缺陷，重構(gòu)的CT方程和原CT方程求解的結(jié)果差別不大，無(wú)缺陷處與理論值相差較小，有缺陷處與理論值相差較大，但是都能通過(guò)概率法準(zhǔn)確地判斷出缺陷的位置。

3.2 算例2

如圖12所示，設(shè)計(jì)一個(gè)數(shù)值模型截面尺寸為800 mm×800 mm，V缺陷=3 500 m/s的兩缺陷數(shù)值模型，其中缺陷尺寸分別為200 mm×200 mm、100 mm×100 mm，將其劃分為8×8的方形網(wǎng)格。通過(guò)計(jì)算均值、標(biāo)準(zhǔn)差、相對(duì)誤差進(jìn)行對(duì)比，其對(duì)比結(jié)果如表3所示。

圖12 兩缺陷數(shù)值模型

方程類(lèi)別均值/(m·s-1)有缺陷無(wú)缺陷標(biāo)準(zhǔn)差/(m2·s-1)有缺陷無(wú)缺陷相對(duì)誤差/%有缺陷無(wú)缺陷彈性波CT方程387244500167．010．643．22彈性波CT方程+8線(xiàn)性變換方程3812447112．8159．48．922．85彈性波CT方程+32線(xiàn)性變換方程3809447210．5156．68．832．84

表2 100 mm×100 mm缺陷計(jì)算結(jié)果對(duì)比

表3 計(jì)算結(jié)果對(duì)比

經(jīng)對(duì)比發(fā)現(xiàn)：重構(gòu)的CT方程和原CT方程求解的結(jié)果差別不大，無(wú)缺陷處與理論值相差較小，有缺陷處與理論值相差較大，但是都能通過(guò)概率法準(zhǔn)確地判斷出缺陷的位置。

4 混凝土試件試驗(yàn)

利用如圖13所示的混凝土試件進(jìn)行實(shí)驗(yàn)，試件尺寸為800 mm×800 mm，其中中心部分為缺陷，大小為200 mm×200 mm。將該測(cè)區(qū)離散化為100 mm×100 mm的網(wǎng)格單元，分別在長(zhǎng)度方向的兩側(cè)布置發(fā)射/接收換能器，相鄰兩個(gè)發(fā)射/接收換能器的間距為100 mm，獲得8×8個(gè)走時(shí)。通過(guò)原始方法和本文方法得到的試件波速數(shù)據(jù)分別如表4、5所示。

圖13 混凝土試件

(m·s-1)

表5 試件A2本文方法彈性波波速 (m·s-1)

基于彈性波CT走時(shí)方程的層析成像結(jié)果以及相關(guān)性重構(gòu)方程的層析成像結(jié)果分別如圖14(a)(b)所示。

圖14 彈性波CT方程重構(gòu)層析成像結(jié)果

通過(guò)兩種方法對(duì)待測(cè)混凝土試件檢測(cè)得到的層析成像進(jìn)行對(duì)比可以發(fā)現(xiàn)：利用重構(gòu)方程得到的圖像對(duì)比度更高，缺陷識(shí)別情況優(yōu)于傳統(tǒng)的理論走時(shí)方程。從而證明了本文算法的穩(wěn)定性更好，同時(shí)對(duì)于偽像的抑制效果更為明顯。

5 結(jié)束語(yǔ)

利用一維傳感器陣列相關(guān)性，提出了一種對(duì)CT層析成像走時(shí)方程的重構(gòu)方法，將線(xiàn)性靜定方程轉(zhuǎn)變?yōu)槌o定方程進(jìn)行優(yōu)化求解；通過(guò)數(shù)值模擬與實(shí)驗(yàn)驗(yàn)證了該方法具有更好的收斂速度和精度，在一定程度上抑制了偽像的產(chǎn)生，從而使反演成像的效果更好。

[1] 徐鈺,段衛(wèi)星,徐維秀,等.高精度初至自動(dòng)拾取綜合方法研究[J].物探與化探,2010(5):595-599.

[2] 左國(guó)平,王彥春,隋榮亮.利用能量比法拾取地震初至的一種改進(jìn)方法[J].石油物探,2004(4):345-347+5.

[3] 向陽(yáng).基于互相關(guān)延時(shí)估計(jì)的波速估計(jì)方法[J].武漢理工大學(xué)學(xué)報(bào)(信息與管理工程版),2003(5):63-65.

[4] 張軍華,趙勇,趙愛(ài)國(guó),等.用小波變換與能量比方法聯(lián)合拾取初至波[J].物探化探計(jì)算技術(shù),2002(4):309-312，336.

[5] 王春艷,樊官民,孟杰.基于廣義互相關(guān)函數(shù)的聲波陣列時(shí)延估計(jì)算法[J].電聲技術(shù),2010(8):37-39.

[6] COPPENS F.First arrival picking on common-offset trace collections for automatic estimation of static correction[J].Geophysical Prospecting,1985,33(8):1212-1232.

[7] 唐小明,吳昊,劉志坤.基于廣義互相關(guān)算法的時(shí)延估計(jì)研究[J].電聲技術(shù),2009(8):71-74.

[8] SUN Hongmei,JIA Ruisheng,DU Qianqian,et al.Cross-correlation analysis and time delay estimation of a homologous micro-seismic signal based on the Hilbert-Huang transform[J].Computers and Geosciences,2016,91(C):98-104.

[9] 王珂,肖鵬峰,馮學(xué)智,等.基于改進(jìn)二維離散希爾伯特變換的圖像邊緣檢測(cè)方法[J].測(cè)繪學(xué)報(bào),2012(3):421-427，433.

[10] 向陽(yáng),彭勇.基于希爾伯特變換的應(yīng)力波波速估計(jì)的研究[J].公路交通科技,2004(1):54-57.

[11] FENG Wang,MEI Quanliu,JIANG Weifan.The performance correlation hilbert time-delay estimation for passive detection[J].Applied Mechanics and Materials,2014,602/605:1768-1771.

MethodofCTImageReconstructionBasedonCorrelationofReceiverArray

ZHANG Yu, CHENG Hua, HUANG Xiaohan, YUAN Ye

(Department of Military Civil Engineering, Logistical Engineering University， Chongqing 401331， China)

Elastic wave CT tomographic imaging technique can extract information of internal defects without damaging the internal structure of concrete members, and reconstruction image by some certain algorithm. This paper proposed a new algorithm of reconstruction equation which took the correlation of the receiver array into consideration, composite equation through linear transformation. Therefore, it turns the linear static equation into statically indeterminate equation and reduces the condition number of coefficient matrix, in order to achieve the objectives of reducing the measurement error and suppressing noise. The numerical simulation shows this method can obtain better images in the same condition, and the experimental evidence suggests this method has low noise and high accuracy in the process of defect discrimination.

elastic wave; inversion algorithm; refactor equation; image reconstruction

2017-06-01

張彧(1993—),女,遼寧人,碩士研究生,主要從事軍事特種結(jié)構(gòu)及檢測(cè)加固研究，E-mail：yu_jade21@163.com；通訊作者 程華(1958—),男,浙江人,教授，主要從事軍事特種結(jié)構(gòu)及檢測(cè)加固研究，E-mail：chwjct@163.com。

張彧，程華，黃曉寒，等.基于陣列傳感器相關(guān)性的CT圖像重建研究[J].重慶理工大學(xué)學(xué)報(bào)(自然科學(xué))，2017(12):181-188.

formatZHANG Yu, CHENG Hua, HUANG Xiaohan, et al.Method of CT Image Reconstruction Based on Correlation of Receiver Array[J].Journal of Chongqing University of Technology(Natural Science)，2017(12):181-188.

10.3969/j.issn.1674-8425(z).2017.12.031

TU391.4;O242

A

1674-8425(2017)12-0181-08

(責(zé)任編輯楊黎麗)