基于陣列傳感器相關性的CT圖像重建研究

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

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

基于陣列傳感器相關性的CT圖像重建研究

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

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

彈性波CT層析成像技術可以在不損傷混凝土構件內部結構的情況下對其內部缺陷的信息進行提取，再利用相關算法重建圖像。對接收點陣列間的相關性加以考慮，通過線性變換的方法組合方程，進而得到了一種重構方程的反演算法。由此將線性靜定方程轉變為超靜定方程，減小系數矩陣的條件數，從而達到了減小測量誤差抑制噪聲影響的目的。數值模擬結果表明:在相同條件下，反演得到的圖像更為清晰。從實驗數據上看，重建所得的圖像噪聲更小，對缺陷判別的準確度更高。

彈性波；反演算法；方程重構；圖像重建

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

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

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

圖1 一維傳感器陣列相關性示意圖

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

圖2 信號相關性示意圖

1 多時窗能量比法

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

(1)

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

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

圖3 滑動時窗能量比法示意圖

在實際工程應用中，信號的采集具有固定采樣頻率，所以對式(1)做離散化處理后得到

(2)

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

(3)

式中：A為該信號道的相對能量；ω為該信號道的樣點總數；α為穩定系數，通常取0.5～2.0。

滑動時窗能量比法僅僅是通過尋找能量比的最大值以獲得初至走時，這種方法存在固有缺陷，因為初至波處的能量比不一定達到最大，從而導致初至走時的拾取錯誤。因此，本文采用多時窗能量比法時不僅考慮了能量比極大值，同時也考慮了能量比次級值，對于低信噪比信號初至走時的拾取精度較好，穩定性也較好。

多時窗能量比法采用4個滑動時窗，如圖4所示，其中L1、L2、L3、L4分別為第1～4時窗。

第1時窗與第2時窗能量比為

(4)

第1時窗與第3時窗能量比為

(5)

第4時窗與第2時窗能量比為

(6)

式中：t為當前采樣點；d為后時窗與第3時窗的間隔；α為穩定系數，一般取值為0.5～2.0。

圖4 多時窗能量比法示意圖

計算步驟為：

1) 根據式(4)計算A1(t)，搜索A1(t)極大值Amax，并記錄其對應的時間T1。

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

圖5 彈性波響應信號

圖6 多時窗能量比法拾取起跳點

圖7 處理后信號

2 廣義二次互相關的希爾伯特變換

2.1 二次互相關

假設信號模型為：

(7)

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

x1(n)和x2(n)的相關函數表示為

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)

假設信號與噪聲以及噪聲與噪聲之間都不相關，則式(8)變為

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

(9)

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

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)

彈性波信號與噪聲的相關函數可以忽略，則式(8)變為

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

(11)

式中：RRS為純信號做二次相關；RRN為噪聲做二次自相關。對于式(9)，當τ=D時，RRR(τ)取得最大值，其峰值點對應的橫坐標為時間延遲。二次相關法的優勢在于能夠進一步降低噪聲對信號的影響，可以在更低信噪比環境下估計時間延遲。

2.2 廣義互相關

廣義互相關算法[7-8]是利用兩個信號的廣義相關函數來估計時延。為了獲得更好的時延估計精度，先利用頻域加權函數對信號進行預濾波，然后再進行相關性的計算，檢測出相關函數峰值的橫坐標即為時延。定義經過加權的廣義相關函數為

(12)

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

(13)

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

圖8 廣義二次互相關算法流程

2.3 希爾伯特變換

希爾波特變換定義如下

(14)

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

圖9 希爾伯特變換示意圖

2.4 方程重構的實現與求解

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

ΔDRm-RnS=ΔTmaxRm-Rn

(15)

式中：ΔDRm-Rn表示相同激勵下第m個接收點和第n個接收點的系數差；S表示慢度向量； ΔTmaxRm-Rn表示相同激勵下第m個接收點和第n個接收點相關函數的極大值。

對添加方程的數量進行了討論：實測過程當中，方程數量的增加會大大增加計算量并帶來誤差，進而影響結果的精度，所以對于8×8的網格單元，附加8個方程進行分析，補充后的方程表達式如下：

圖10 一維傳感器陣列相關性示意圖

3 數值算例

3.1 算例1

如圖11所示，設計2個模型截面尺寸為800 mm×800 mm，缺陷尺寸為200 mm×200 mm、100 mm×100 mm的單缺陷數值模型，并將其劃分為8×8的方形網格。通過計算均值、標準差、相對誤差進行對比，其對比結果分別如表1、2所示。

圖11 不同大小缺陷數值模型

經對比發現：對于不同尺寸的缺陷，重構的CT方程和原CT方程求解的結果差別不大，無缺陷處與理論值相差較小，有缺陷處與理論值相差較大，但是都能通過概率法準確地判斷出缺陷的位置。

3.2 算例2

如圖12所示，設計一個數值模型截面尺寸為800 mm×800 mm，V缺陷=3 500 m/s的兩缺陷數值模型，其中缺陷尺寸分別為200 mm×200 mm、100 mm×100 mm，將其劃分為8×8的方形網格。通過計算均值、標準差、相對誤差進行對比，其對比結果如表3所示。

圖12 兩缺陷數值模型

方程類別均值/(m·s-1)有缺陷無缺陷標準差/(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缺陷計算結果對比

表3 計算結果對比

經對比發現：重構的CT方程和原CT方程求解的結果差別不大，無缺陷處與理論值相差較小，有缺陷處與理論值相差較大，但是都能通過概率法準確地判斷出缺陷的位置。

4 混凝土試件試驗

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

圖13 混凝土試件

(m·s-1)

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

基于彈性波CT走時方程的層析成像結果以及相關性重構方程的層析成像結果分別如圖14(a)(b)所示。

圖14 彈性波CT方程重構層析成像結果

通過兩種方法對待測混凝土試件檢測得到的層析成像進行對比可以發現：利用重構方程得到的圖像對比度更高，缺陷識別情況優于傳統的理論走時方程。從而證明了本文算法的穩定性更好，同時對于偽像的抑制效果更為明顯。

5 結束語

利用一維傳感器陣列相關性，提出了一種對CT層析成像走時方程的重構方法，將線性靜定方程轉變為超靜定方程進行優化求解；通過數值模擬與實驗驗證了該方法具有更好的收斂速度和精度，在一定程度上抑制了偽像的產生，從而使反演成像的效果更好。

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

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

[3] 向陽.基于互相關延時估計的波速估計方法[J].武漢理工大學學報(信息與管理工程版),2003(5):63-65.

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

[5] 王春艷,樊官民,孟杰.基于廣義互相關函數的聲波陣列時延估計算法[J].電聲技術,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] 唐小明,吳昊,劉志坤.基于廣義互相關算法的時延估計研究[J].電聲技術,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] 王珂,肖鵬峰,馮學智,等.基于改進二維離散希爾伯特變換的圖像邊緣檢測方法[J].測繪學報,2012(3):421-427，433.

[10] 向陽,彭勇.基于希爾伯特變換的應力波波速估計的研究[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—),女,遼寧人,碩士研究生,主要從事軍事特種結構及檢測加固研究，E-mail：yu_jade21@163.com；通訊作者 程華(1958—),男,浙江人,教授，主要從事軍事特種結構及檢測加固研究，E-mail：chwjct@163.com。

張彧，程華，黃曉寒，等.基于陣列傳感器相關性的CT圖像重建研究[J].重慶理工大學學報(自然科學)，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

(責任編輯楊黎麗)