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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

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

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

在實(shí)際工程應(yīng)用中，信號的采集具有固定采樣頻率，所以對式(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ā)生，對式(2)分子分母同時(shí)加上穩(wěn)定因子，得

(3)

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

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

多時(shí)窗能量比法采用4個(gè)滑動時(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，并記錄其對應(yīng)的時(shí)間T1。

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

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

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

圖7 處理后信號

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

2.1 二次互相關(guān)

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

(7)

其中：s(n)為彈性波信號；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è)信號與噪聲以及噪聲與噪聲之間都不相關(guān)，則式(8)變?yōu)?/p>

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

(9)

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

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

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

(11)

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

2.2 廣義互相關(guān)

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

(12)

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

(13)

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

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

2.3 希爾伯特變換

希爾波特變換定義如下

(14)

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

圖9 希爾伯特變換示意圖

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

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

ΔDRm-RnS=ΔTmaxRm-Rn

(15)

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

對添加方程的數(shù)量進(jìn)行了討論：實(shí)測過程當(dāng)中，方程數(shù)量的增加會大大增加計(jì)算量并帶來誤差，進(jìn)而影響結(jié)果的精度，所以對于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)格。通過計(jì)算均值、標(biāo)準(zhǔn)差、相對誤差進(jìn)行對比，其對比結(jié)果分別如表1、2所示。

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

經(jīng)對比發(fā)現(xiàn)：對于不同尺寸的缺陷，重構(gòu)的CT方程和原CT方程求解的結(jié)果差別不大，無缺陷處與理論值相差較小，有缺陷處與理論值相差較大，但是都能通過概率法準(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)格。通過計(jì)算均值、標(biāo)準(zhǔn)差、相對誤差進(jìn)行對比，其對比結(jié)果如表3所示。

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

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

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

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

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

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

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

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

5 結(jié)束語

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

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

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

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

[4] 張軍華,趙勇,趙愛國,等.用小波變換與能量比方法聯(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)二維離散希爾伯特變換的圖像邊緣檢測方法[J].測繪學(xué)報(bào),2012(3):421-427，433.

[10] 向陽,彭勇.基于希爾伯特變換的應(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)及檢測加固研究，E-mail：yu_jade21@163.com；通訊作者 程華(1958—),男,浙江人,教授，主要從事軍事特種結(jié)構(gòu)及檢測加固研究，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é)任編輯楊黎麗)