基于稀疏重構(gòu)的尾波干涉成像方法*

張濤 侯宏? 鮑明

1) (西北工業(yè)大學(xué)航海學(xué)院, 海洋聲學(xué)信息感知工業(yè)和信息化部重點(diǎn)實(shí)驗(yàn)室, 西安 710072)

2) (中國(guó)科學(xué)院噪聲與振動(dòng)重點(diǎn)實(shí)驗(yàn)室(聲學(xué)研究所), 北京 100190)

尾波干涉成像是利用尾波時(shí)延和擴(kuò)散近似敏感核來反演散射介質(zhì)中微小速度擾動(dòng)空間分布的技術(shù).該問題本質(zhì)上是一個(gè)欠定問題, 一般無確定解, 導(dǎo)致其難以精確定位介質(zhì)中微小波速變化的區(qū)域.為解決上述缺陷, 本文利用速度擾動(dòng)分布在空間上具有稀疏性的特點(diǎn), 提出了一種基于壓縮感知理論中稀疏重構(gòu)算法的尾波干涉成像方法.該方法可在散射介質(zhì)中較準(zhǔn)確地獲取速度擾動(dòng)的空間位置和范圍, 同時(shí)具有較高的計(jì)算效率.數(shù)值仿真和實(shí)驗(yàn)結(jié)果表明:相比于現(xiàn)有的線性最小二乘差分成像方法, 無論是單個(gè)還是多個(gè)擾動(dòng)區(qū)域,該方法均能更精確地進(jìn)行定位成像, 同時(shí)明顯減少了計(jì)算時(shí)間.

1 引 言

尾波干涉 (coda wave interferometry, CWI),是一種利用介質(zhì)中的多次散射波檢測(cè)介質(zhì)微小變化的測(cè)量方法.多次散射波是聲信號(hào)或振動(dòng)信號(hào)在散射介質(zhì)中傳播時(shí), 由于介質(zhì)的非均勻性(存在如顆粒、裂縫、包體以及孔洞等非均質(zhì)體)而發(fā)生多次散射、反射產(chǎn)生的, 在波形圖中顯示為直達(dá)波后面拖著的長(zhǎng)長(zhǎng)的“尾巴”, 因此也稱尾波.尾波由于在介質(zhì)中來回傳播, 對(duì)介質(zhì)更密集地采樣, 使得介質(zhì)中的微小變化不斷迭加放大, 從而對(duì)介質(zhì)性質(zhì)的微弱改變具有極高的敏感性, 檢測(cè)出直達(dá)波無法識(shí)別的微小變化.自2002年CWI法首次由Snieder等[1?3]提出后, 便在地球物理學(xué)和材料科學(xué)等領(lǐng)域得到了廣泛的研究與應(yīng)用[4?9].文獻(xiàn)[7]在SAFOD(San Andreas Fault Observatory at Depth) 進(jìn)行了一次主動(dòng)震源井間實(shí)驗(yàn), 觀測(cè)到距試驗(yàn)點(diǎn)3 km外的一個(gè)3級(jí)地震發(fā)生前數(shù)小時(shí), 井內(nèi)波速增加了0.8%.文獻(xiàn)[8]運(yùn)用CWI測(cè)得水泥樣品水飽和進(jìn)程和樓房承載柱的尾波波速隨環(huán)境濕度的變化,水飽和使得尾波波速相對(duì)變化0.7%, 樓房承載柱隨相對(duì)濕度每變化1%, 尾波波速變化0.213%.文獻(xiàn)[9]基于CWI的測(cè)量方法, 利用超聲尾波在實(shí)驗(yàn)室觀測(cè)了1.5 m巖石斷層的黏滑過程, 獲得了精度高達(dá) 10–6的相對(duì)波速變化, 相當(dāng)于~10 kPa 的應(yīng)力變化.以上研究都取得了一定成果, 均通過CWI的測(cè)量方法獲取了介質(zhì)的相對(duì)波速變化量, 但是這些變化量?jī)H僅是觀測(cè)介質(zhì)各位置相對(duì)波速變化的加權(quán)平均, 并不能反映波速變化的空間分布情況.然而, 正確定位波速變化區(qū)域?qū)τ诜治鲆鸩ㄋ僮兓奈锢頇C(jī)制十分重要, 因此對(duì)波速變化的空間分布成像近年來成為了CWI法的研究前沿與難點(diǎn)之一.早在 2005年, Pacheco 和 Snieder[10]便 對(duì)CWI測(cè)量方法進(jìn)行了延伸, 利用多次散射波場(chǎng)的擴(kuò)散近似, 提出了一種通過擴(kuò)散近似敏感核將平均走時(shí)延遲與介質(zhì)波速局部變化聯(lián)系起來的技術(shù), 使定位局部微小變化成為了可能.文獻(xiàn)[11]分析了2010年6月至 12月La Reunion 島上 Piton de la Fournaise活火山的連續(xù)環(huán)境噪聲記錄, 利用10月的大爆發(fā)前觀察到的速度變化確定了爆發(fā)點(diǎn)速度顯著下降(高達(dá)0.6%)的區(qū)域.文獻(xiàn)[12]使用環(huán)境噪聲互相關(guān)和延展法處理了墨西哥Volcán de Colima地區(qū)超過15年的連續(xù)記錄地震數(shù)據(jù), 最明顯的速度變化是2003年Tecomán附近的7.4級(jí)地震期間下降了2.6%, 并采用了基于輻射傳輸近似的敏感核反演定位水平面上的擾動(dòng).文獻(xiàn)[13]描述了一種LOCADIFF技術(shù), 利用CWI法以及結(jié)合擴(kuò)散傳播模型的最大似然法反演定位非均勻強(qiáng)散射介質(zhì)中的微小變化, 并通過二維有限差分進(jìn)行了數(shù)值模擬.文獻(xiàn)[14]將LOCADIFF這種超聲波成像技術(shù)成功應(yīng)用在混凝土結(jié)構(gòu)無損檢測(cè)與評(píng)價(jià)中.對(duì)于CWI結(jié)合敏感核成像介質(zhì)波速變化的空間分布以定位局部微小改變的研究, 國(guó)內(nèi)外目前涉及較少, 反演問題是這種成像技術(shù)發(fā)展與應(yīng)用的阻力之一.尾波干涉成像技術(shù)歸根到底是對(duì)一個(gè)矩陣形式的反問題進(jìn)行求解, 該問題是一個(gè)欠定問題, 有無窮多解, 同時(shí)具有病態(tài)性, 難以精確求解.前述研究中均采用地球物理反演中常用的線性最小二乘差分反演方法, 該方法存在求解精度不高、參數(shù)選取困難的缺陷, 并且難以定位多個(gè)擾動(dòng)點(diǎn), 在實(shí)際工程中的應(yīng)用十分困難.

另一方面, 由 Donoho[15]以 及 Candes 和Romberg[16]于2006年正式提出的壓縮感知理論(compressed sensing, CS)近年來逐漸興起并趨于成熟, 在雷達(dá)探測(cè)、圖像處理、聲吶定位等領(lǐng)域都有了廣泛的應(yīng)用, 在物理成像方面也表現(xiàn)出強(qiáng)大的生命力[17?21].文獻(xiàn) [17]采用分塊壓縮感知思想, 提出一種基于lp范數(shù)的壓縮感知圖像重建算法, 提高重建圖像質(zhì)量的同時(shí)縮短了成像時(shí)間.文獻(xiàn)[18]提出了基于稀疏測(cè)量的超分辨壓縮感知鬼成像重建模型, 并搭建了實(shí)驗(yàn)裝置, 驗(yàn)證了該模型的有效性與優(yōu)越性, 推動(dòng)了鬼成像的實(shí)用進(jìn)展.文獻(xiàn)[19]利用聲源的空間稀疏性, 在壓縮感知理論下建立了矢量陣聚焦定位新方法, 克服了相干聲源分辨困難、貢獻(xiàn)評(píng)價(jià)不準(zhǔn)確等問題, 且定位精度高, 背景抑制能力強(qiáng).CS理論指出, 若某信號(hào)是可壓縮的或經(jīng)某種變換后稀疏, 則可通過隨機(jī)觀測(cè)矩陣(與變換基不相關(guān))將高維信號(hào)投影至低維空間中, 最終通過求解優(yōu)化問題從少量低維數(shù)據(jù)高概率地重構(gòu)原始高維信號(hào).而在尾波干涉成像問題中, 實(shí)際擾動(dòng)通常是空間上聚焦的, 擾動(dòng)的空間區(qū)域相比于監(jiān)測(cè)的整體介質(zhì)區(qū)域往往是稀疏的, 恰好滿足CS理論對(duì)于信號(hào)稀疏性的要求.本文基于CS理論, 提出了一種尾波干涉成像的新的解算方法.該方法首先通過CWI和擴(kuò)散波敏感核構(gòu)建了用于成像的反演模型; 其次引入CS理論框架下的稀疏變換將欠定方程重新構(gòu)造; 最終通過壓縮感知重構(gòu)算法——正交匹配追蹤算法 (orthogonal matching pursuit,OMP)[22]求解該欠定方程, 并獲得速度擾動(dòng)的空間成像.此方法與先前的線性最小二乘差分成像方法相比, 在單擾動(dòng)和多擾動(dòng)情況下對(duì)于擾動(dòng)區(qū)域位置和范圍的確定都更加準(zhǔn)確, 且執(zhí)行簡(jiǎn)易, 沒有復(fù)雜的參數(shù)選取, 同時(shí)還具有更高的計(jì)算效率.最后結(jié)合具體仿真算例, 給出了與線性最小二乘差分成像方法的比較結(jié)果, 證明了該成像方法的可行性、有效性和高效率.

2 尾波干涉成像基本原理

2.1 尾波干涉基本原理

尾波干涉理論指出, 多散射介質(zhì)中的微小改變將導(dǎo)致波速的微弱變化, 主要表現(xiàn)為尾波信號(hào)在擾動(dòng)前后的相位偏移.圖1顯示了微小結(jié)構(gòu)變化發(fā)生前后的多散射介質(zhì)中某處記錄的波形信號(hào), 插圖顯示了直達(dá)波(圖1(a))和尾波(圖1(b))波形的詳細(xì)信息.不難看出, 直達(dá)波于擾動(dòng)前后在幅值和相位上幾乎沒有變化, 而尾波則表現(xiàn)出明顯的走時(shí)延遲, 該時(shí)延即為尾波對(duì)于介質(zhì)中微弱變化多次采樣后的體現(xiàn), 是直達(dá)波所不能包含的物理信息.

尾波干涉理論的實(shí)質(zhì)是通過測(cè)量?jī)蓚€(gè)不同時(shí)刻尾波列的相位差來獲取介質(zhì)在該時(shí)段的波速變化[1?3].互相關(guān)法和延展法是目前基于此理論來提取介質(zhì)波速變化的主要方法[23,24], 考慮到仿真信號(hào)信噪比較高, 本文采用計(jì)算效率更高的互相關(guān)法[23].步驟如下:假設(shè)已獲取介質(zhì)變化前后的兩列尾波信號(hào)uunp和uper, 通過兩列波的互相關(guān)計(jì)算獲得波速變化, 即

圖1 多散射介質(zhì)中擾動(dòng)前后波形的比較 (a) 直達(dá)波擾動(dòng)前后的波形; (b)尾波擾動(dòng)前后的波形Fig.1.Comparison between typical time traces of a wave propagating in a multiple scattering medium before and after a small perturbation:(a) The first arrival waves before and after a small perturbation; (b) the coda waves before and after a small perturbation.

式中, 2T表示要分析的尾波部分的時(shí)間窗長(zhǎng)度;t為時(shí)間窗中心位置;ts表示互相關(guān)中所用的走時(shí)差;R的值表示兩列波的相似程度.

當(dāng)ts的某取值tsmax使得R(ts) 取最大值時(shí), 該走時(shí)差tsmax即為所分析尾波部分?jǐn)_動(dòng)前后的時(shí)間延遲δt.那么根據(jù)Snieder[1]對(duì)尾波干涉的系統(tǒng)闡述, 可得兩列信號(hào)的相對(duì)波速變化

式中 δυ表示波速擾動(dòng), 是整體介質(zhì)波速變化的加權(quán)平均;υ為擾動(dòng)前的介質(zhì)波速.

2.2 敏感核

為了獲取局部擾動(dòng)的空間分布情況, 需要將上述尾波干涉求得的總體平均時(shí)延δt與介質(zhì)波速局部變化聯(lián)系起來.然而, 由于介質(zhì)的復(fù)雜性, 多散射效應(yīng)產(chǎn)生的尾波的傳播路徑長(zhǎng)而無序, 因此難以將每個(gè)時(shí)延與一組特定的傳播軌跡配對(duì), 確定變化的具體來源也就十分困難.一些學(xué)者[10,13]并沒有嘗試確定每個(gè)尾波序列的傳播路徑, 而是研究了一種基于時(shí)空傳播統(tǒng)計(jì)描述的替代方法, 引入敏感核,K(s1,s2,x0,t) 來描述波列在位置s1發(fā)射, 途經(jīng)位置x0, 并在總傳播時(shí)間t后在位置s2被接收的概率大小, 此概率描述了尾波在位置x0處傳播時(shí)間的空間密度.

式中的p(s1,s2,t) 表示波列從s1經(jīng)過時(shí)間t到達(dá)s2的概率.在實(shí)際應(yīng)用中, 此概率可以用波場(chǎng)強(qiáng)度等效, 波場(chǎng)強(qiáng)度是位置和時(shí)間的函數(shù), 可以通過擴(kuò)散方程的解來近似.在無限大二維介質(zhì)中, 擴(kuò)散方程的解為

式中, |s1?s2| 表示s1和s2之間的距離;D為介質(zhì)的散射系數(shù), 與其物理性質(zhì)有關(guān).通過激勵(lì)源與接收點(diǎn)的位置和尾波的傳播時(shí)間, 對(duì)介質(zhì)空間逐點(diǎn)計(jì)算K(s1,s2,x0,t) , 即可得到敏感核的空間分布.圖2所示為尾波觀察時(shí)間t分別為1 s和5 s, 散射系數(shù)為 5.8×105m2/s , 激勵(lì)源 (S)與接收點(diǎn) (R)間距5000 m時(shí)的敏感核空間分布圖.從圖2中可以看出, 敏感核在三維空間上呈馬鞍形分布, 在激勵(lì)點(diǎn)和接收點(diǎn)處達(dá)到峰值, 敏感核的值隨著與這兩點(diǎn)的距離增大而降低, 敏感核的空間分布也隨著尾波觀察時(shí)間的增加而展寬.t= 1 s時(shí), 以直達(dá)波和單散射波為主, 對(duì)介質(zhì)的采樣密度和范圍較小,t= 5 s時(shí), 以多次散射產(chǎn)生的尾波為主, 對(duì)介質(zhì)進(jìn)行密集采樣的同時(shí)增加了探測(cè)范圍.

2.3 反演成像模型

基于擴(kuò)散近似敏感核, 文獻(xiàn)[10,11,13]提出了一個(gè)線性模型, 將尾波干涉法獲取的平均時(shí)延數(shù)據(jù)通過敏感核項(xiàng)K(s1,s2,x0,t) 與相應(yīng)的局部變化聯(lián)系起來.

圖2 基于擴(kuò)散近似的二維敏感核示例 (a) t= 1 s時(shí)的敏感核空間分布; (b) t= 5 s 時(shí)的敏感核空間分布; (c) t= 1 s 時(shí)的敏感核俯視圖; (d) t= 5 s時(shí)的敏感核俯視圖Fig.2.Examples of sensitivity kernel based on the diffusion approximation in 2-D:(a) Spatial representation of the sensitivity kernel when t= 1 s; (b) spatial representation of the sensitivity kernel when t= 5 s; (c) vertical view of the sensitivity kernel when t=1 s; (d) vertical view of the sensitivity kernel when t= 5 s.

為了使探測(cè)范圍盡量覆蓋整個(gè)介質(zhì), 尾波干涉成像需要大量的敏感核, 因此會(huì)有許多形如(5)式的方程, 為清晰表達(dá), 可將線性方程組寫成矩陣形式:

式中,T∈RM,1表示通過不同的激勵(lì)接收對(duì)和不同的尾波觀察窗口獲取的M個(gè)時(shí)延數(shù)據(jù);V∈RN,1表示空間各網(wǎng)格點(diǎn)對(duì)應(yīng)的相對(duì)波速擾動(dòng)值,N為網(wǎng)格總數(shù);G∈RM,N表示T到V的映射,其矩陣元素Gij=?sKij, 其中Kij為敏感核矩陣K的元素,K∈RM,N的每行代表一個(gè)敏感核, 每列代表一個(gè)網(wǎng)格, 每個(gè)網(wǎng)格的面積為 ?s.

利用線性最小二乘差分法對(duì)V值進(jìn)行反演, 即

式中,V0為先驗(yàn)初始值, 一般為零矩陣;GT表示G的轉(zhuǎn)置;Cd為測(cè)得數(shù)據(jù)的對(duì)角協(xié)方差矩陣;Cm為介質(zhì)空間元素平滑矩陣, 可用下式計(jì)算:

式中,σm為先驗(yàn)標(biāo)準(zhǔn)差;λ0為網(wǎng)格間距;λ為相關(guān)長(zhǎng)度.σm和λ一般通過 L 曲線法[25]選取.

在反演計(jì)算時(shí), 還需要循環(huán)迭代Cm以尋找最優(yōu)的V值, 迭代公式如下:

3 基于稀疏重構(gòu)的尾波干涉成像方法

壓縮感知理論指出, 假定N維真實(shí)信號(hào)x∈RN在某變換域下是L稀疏的(即則可利用信號(hào)x的任意M(M≈Lln(N/L)) 個(gè)線性測(cè)量值高概率精確重構(gòu)x.壓縮感知原理主要包括三個(gè)過程:信號(hào)的稀疏表示、測(cè)量矩陣構(gòu)建以及算法重構(gòu).

首先, 對(duì)信號(hào)進(jìn)行稀疏表示如下:

式中,ψ∈RN為稀疏變換矩陣;a∈RN,1為x的稀疏表示.

其次, 構(gòu)建一個(gè)與變換基底ψ不相關(guān)的觀測(cè)矩陣對(duì)a進(jìn)行線性觀測(cè)

式中,?∈RM,N為觀測(cè)矩陣;y∈RM,1是a在觀測(cè)矩陣下的觀測(cè)值;Θ=?ψ稱為傳感矩陣.

最后, 對(duì)信號(hào)進(jìn)行重構(gòu), 是壓縮感知理論中極為關(guān)鍵的一環(huán).(11) 式中y=Θx稱為欠定方程,y中的未知量有N個(gè), 而方程僅有M個(gè), 理論上該反問題有無窮多個(gè)解, 但如果Θ滿足有限等距性質(zhì) (restricted isometry property, RIP), 則有唯一最稀疏解.該過程等效為求解一個(gè)最優(yōu)化問題

求解 (12)式是最小l0范數(shù)優(yōu)化問題, 屬于NP-hard問題, 直接求解此問題數(shù)值極不穩(wěn)定, 且計(jì)算復(fù)雜度高, 實(shí)際工程中實(shí)施困難.Candes[26]提出,在一定條件下, 可以利用l1范數(shù)代替l0范數(shù), 即

因此, (12)式的問題轉(zhuǎn)化為(13)式的凸優(yōu)化問題,可以通過線性規(guī)劃理論求解, 本文采用應(yīng)用廣泛的OMP法作為重構(gòu)算法.

可以觀察到(6)式和(11)式在數(shù)學(xué)形式上相似, 且擾動(dòng)的空間區(qū)域相比于整體介質(zhì)區(qū)域是稀疏的, 恰好滿足壓縮感知理論對(duì)于信號(hào)稀疏性的要求, 因此, 如果將G,T和V分別視為觀測(cè)矩陣、測(cè)量值矩陣和真實(shí)信號(hào), 則可在壓縮感知理論框架下, 通過信號(hào)稀疏表示和稀疏信號(hào)重構(gòu)的方法求解(6)式.

本文引入稀疏變換矩陣P對(duì)V進(jìn)行離散稀疏變換.從而有

式中,GCS=GP?1,VCS=PV.常用的稀疏變換矩陣有離散傅里葉變換矩陣、離散余弦變換矩陣、離散小波變換矩陣等, 圖像信號(hào)在二維小波變換域是稀疏的, 聲信號(hào)在傅里葉變換域是稀疏的, 本文本質(zhì)上是對(duì)一維聲信號(hào)的處理, 因此采用離散傅里葉變換矩陣.利用OMP等凸優(yōu)化算法求解出方程(14)的稀疏解VCS, 最終通過稀疏逆變換即可求得原始解V, 即

整個(gè)過程無須使用P而僅須用到其逆矩陣P?1,對(duì)于傅里葉變換矩陣,P?1=PT.

4 數(shù)值仿真與結(jié)果

為了驗(yàn)證基于稀疏重構(gòu)算法的成像方法, 現(xiàn)通過數(shù)值仿真進(jìn)行分析.本數(shù)值仿真使用的激勵(lì)源為 15 Hz 主頻的雷克子波, 在 1 0km×10km 的非均勻散射介質(zhì)中心處激發(fā), 持續(xù)時(shí)間為 6π/100s ,時(shí)間采樣間隔 ?t為0.001 s.通過傳播速度的不同模擬介質(zhì)的非均勻散射特性, 將介質(zhì)的速度場(chǎng)劃分為 2 00×200 網(wǎng)格, 產(chǎn)生均值為 5 000m/s , 標(biāo)準(zhǔn)差為500m/s的隨機(jī)數(shù)矩陣作為速度矩陣, 其中最大速度為 7 320.95m/s , 最小速度為 3 045.29m/s , 如圖3所示.接收點(diǎn)的布設(shè)原則是確保均勻覆蓋介質(zhì)區(qū)域, 本仿真中采用的36個(gè)接收點(diǎn)和1個(gè)激勵(lì)源的位置分布如圖4所示.

本數(shù)值仿真運(yùn)用二維聲波方程模擬波的傳播,

式中,p(x,y,t) 為質(zhì)點(diǎn)的位移函數(shù),f(x,y,t) 為激勵(lì)源函數(shù),v(x,y) 為介質(zhì)在 (x,y) 處的速度.采用時(shí)間二階、空間八階的有限差分法對(duì)(16)式進(jìn)行數(shù)值離散, 即有

圖3 二維速度場(chǎng)模型Fig.3.2-D velocity field model.

圖4 激勵(lì)源及接收點(diǎn)布設(shè)Fig.4.Layout of the source and receivers.

式中,nx,ny分別為劃分在x,y方向的網(wǎng)格個(gè)數(shù);nt為時(shí)間間隔個(gè)數(shù); ?h,?t分別為空間、時(shí)間網(wǎng)格上的劃分步長(zhǎng), ?h=50m , ?t=0.001s ;pk(i,j)表示第k次(對(duì)應(yīng)的時(shí)間)迭代時(shí)在 (i,j) 處的位移.

為了模擬相對(duì)無限大的傳播空間, 本文仿真的速度場(chǎng)邊界設(shè)置為吸收邊界, 采用Clayton_Engquist_majda吸收邊界算法, 反射率約為2.5%, 滿足本研究的要求[27].

下面基于以上仿真環(huán)境, 根據(jù)第2和第3節(jié)的原理, 給出針對(duì)五個(gè)速度擾動(dòng)算例的線性最小二乘差分成像方法和本文成像方法的性能比較.其中,速度擾動(dòng)的空間分布通過將介質(zhì)剖分為 2 0×20 的網(wǎng)格來離散表達(dá), 則有400個(gè)待反演的未知量, 網(wǎng)格單元邊長(zhǎng)為500 m; 選擇尾波觀察時(shí)段為1.5—4.7 s, 每段窗口長(zhǎng)度為 0.5 s, 重疊 0.2 s, 每個(gè)接收點(diǎn)可以獲取擾動(dòng)前后各10段尾波數(shù)據(jù), 36對(duì)收發(fā)對(duì)即可計(jì)算得到360個(gè)時(shí)延數(shù)據(jù); 敏感核的構(gòu)造中散射系數(shù)D=8×104m2/s ; 線性最小二乘差分方法中的相對(duì)長(zhǎng)度λ=750m.

對(duì)于單點(diǎn)擾動(dòng), 通過算例1和算例2進(jìn)行了兩種算法的成像對(duì)比, 如圖5和圖6所示.其中, 算例1在速度場(chǎng)的左下方(具體位置見圖5藍(lán)色方框)施加一個(gè)大小為+0.5%, 面積為 1 km×2km 的速度擾動(dòng)區(qū)域; 算例2在圖6右上方藍(lán)色方框區(qū)域添加一個(gè)大小為+0.5%, 面積為 2 km×2km 的速度波動(dòng).

圖5(a)和圖6(a)分別為算例1和算例2利用線性最小二乘差分方法迭代10次的結(jié)果, 通過L曲線法獲取的σm分別為 0.00328和 0.00266,圖5(b)和圖6(b)分別為算例1和算例2在壓縮感知理論框架下利用OMP算法迭代5次和7次的結(jié)果.從算例1的結(jié)果可以看出, 兩種算法均能準(zhǔn)確地突顯出擾動(dòng)的位置, 而在擾動(dòng)范圍的確定方面本文方法更為接近真實(shí)值; 但從圖6中可以清晰地看出, 線性最小二乘方法反演出的算例2擾動(dòng)區(qū)域與真實(shí)設(shè)置相距甚遠(yuǎn), 而本文方法依然能得到較好的反演結(jié)果.可見, 本文方法在單擾動(dòng)情況下較線性最小二乘法有更高的精確性與穩(wěn)定性.

圖7展示了實(shí)測(cè)數(shù)據(jù)的處理結(jié)果.該實(shí)驗(yàn)是一個(gè)小型的模擬實(shí)驗(yàn), 利用激振器發(fā)出前述雷克子波, 作為對(duì) 1.5m×1.2m 的石膏板小樣品的激勵(lì)信號(hào), 通過石膏板上布設(shè)的五個(gè)加速度傳感器采集微小擾動(dòng)發(fā)生前后的波形信號(hào).石膏板的波速與石膏板的含水率有關(guān), 含水率增加, 板中超聲波傳播速度相應(yīng)減小, 因此, 本實(shí)驗(yàn)采用加濕的方式對(duì)樣品介質(zhì)添加擾動(dòng).

將石膏板剖分為 1 0×10 的網(wǎng)格, 則待反演量的個(gè)數(shù)為100, 網(wǎng)格單元長(zhǎng)寬分別為 0.15 m和0.12 m, 在石膏板的中心位置加水以增加中心點(diǎn)處的含水率作為局部擾動(dòng); 尾波分析時(shí)段為0.2—2.8 s,窗口長(zhǎng)度均為 0.6 s, 重疊 0.4 s, 共 11 段尾波數(shù)據(jù),則可獲得 5 ×11=55 個(gè)時(shí)延數(shù)據(jù); 敏感核構(gòu)造中的散射系數(shù)D=140m2/s ; 線性最小二乘法中的相對(duì) 長(zhǎng)度λ=0.1m , 網(wǎng)格 間距λ0取0.134 m.

圖5 算例 1 (a) 線性最小二乘法的反演圖像; (b) 本文方法的反演圖像Fig.5.Case 1:(a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

圖6 算例 2 (a) 線性最小二乘法的反演圖像; (b) 本文方法的反演圖像Fig.6.Case 2:(a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

圖7 實(shí)驗(yàn)數(shù)據(jù)處理結(jié)果 (a) 線性最小二乘法的反演圖像; (b) 本文方法的反演圖像Fig.7.The results of experimental data processing:(a) Inversion image of linear least squares method; (b)inversion image of the method in this paper.

在該模擬實(shí)驗(yàn)中, 僅僅對(duì)板的中心處進(jìn)行了加濕, 擾動(dòng)的空間分布已足夠稀疏, 因此在使用稀疏重構(gòu)方法時(shí), 可以直接求解(6)式的最小l1范數(shù),而無需進(jìn)行稀疏變換, 圖7(b)即為稀疏重構(gòu)法的成像結(jié)果, 圖7(a)為線性最小二乘法的迭代結(jié)果,L曲 線 法 求 得σm為 0.33698.對(duì) 比 7(a)和 圖7(b)可以看出, 兩者對(duì)于擾動(dòng)位置均有不錯(cuò)的定位效果, 但稀疏重構(gòu)法的反演結(jié)果中, 擾動(dòng)更為聚焦,范圍更接近于實(shí)際擾動(dòng)的施加情況.

前面討論了單點(diǎn)擾動(dòng)算例, 并驗(yàn)證了所提方法的優(yōu)越性能.在實(shí)際中, 往往會(huì)出現(xiàn)多個(gè)擾動(dòng)的情況, 因此算例3至5考慮了多個(gè)擾動(dòng)區(qū)域的仿真,圖8—圖10分別給出了兩種算法反演結(jié)果的對(duì)比,圖中的藍(lán)色方框和黃色方框分別表示在該區(qū)域內(nèi)施加+0.5%和–0.5%的速度變化.

圖8 算例 3 (a) 線性最小二乘法的反演圖像; (b) 本文方法的反演圖像Fig.8.Case 3:(a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

圖9 算例 4 (a) 線性最小二乘法的反演圖像; (b) 本文方法的反演圖像Fig.9.Case 4:(a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

圖10 算例 5 (a) 線性最小二乘法的反演圖像; (b) 本文方法的反演圖像Fig.10.Case 5:(a) Inversion image of linear least squares method; (b) inversion image of the method in this paper.

圖8(a), 圖9(a)和圖10(a)分別是算例 3, 4,5用線性最小二乘法迭代10次的結(jié)果, L曲線法獲 取 的σm分 別 為 0.01244, 0.00201和 0.00208,圖8(b), 圖9(b)和圖10(b)分別為算例 3, 4, 5 用本文方法迭代5次、5次和4次的結(jié)果.從圖9和圖10可以直觀地看出, 本文方法對(duì)于多擾動(dòng)區(qū)域的反演質(zhì)量明顯優(yōu)于線性最小二乘法, 線性最小二乘法在算例4和算例5的結(jié)果中, 僅能定位出部分?jǐn)_動(dòng), 且其位置與范圍也有一定程度的偏離真值,而本文方法對(duì)于設(shè)定的多擾動(dòng)區(qū)域, 均能較為準(zhǔn)確地反演出其位置與范圍; 而從算例3的反演圖像來看, 圖8(a)確定的擾動(dòng)區(qū)域與真值相差甚遠(yuǎn), 甚至失去了參考價(jià)值, 相比之下, 圖8(b)的結(jié)果更為接近真值, 有利于正確分析引起波速變化的物理機(jī)制.

表1 線性最小二乘法與本文方法計(jì)算成像時(shí)間對(duì)比Table 1.The comparison of imaging time between linear least squares method and the method in this paper.

此外, 兩種方法在5個(gè)算例下反演成像時(shí)間對(duì)比如表1所列.

從表1可以看出, 在5個(gè)算例中, 本文方法的成像時(shí)間均小于線性最小二乘法, 平均減少了86%左右.同時(shí), 模擬實(shí)驗(yàn)中兩種方法的成像時(shí)間分別為 1.075540 s和 0.122640 s, 稀疏重構(gòu)法的成像時(shí)間比線性最小二乘法少了88.6%.

綜合以上數(shù)值仿真及實(shí)驗(yàn)結(jié)果可知, 壓縮感知理論下的尾波干涉成像方法的性能優(yōu)于線性最小二乘差分法, 不僅在多擾動(dòng)識(shí)別能力、定位精度上有所提高, 成像時(shí)間也大為縮短.

5 結(jié) 論

在尾波干涉成像中, 現(xiàn)有的線性最小二乘差分成像方法定位精度低, 且在多擾動(dòng)區(qū)域的識(shí)別中幾乎失效.針對(duì)該問題, 充分利用擾動(dòng)區(qū)域的空間稀疏性, 提出了一種基于稀疏重構(gòu)的尾波干涉成像新方法.該方法首先根據(jù)尾波干涉法獲取的時(shí)延數(shù)據(jù)以及擴(kuò)散近似敏感核矩陣建立反演成像的欠定方程, 其次在壓縮感知框架下, 通過信號(hào)的稀疏表示進(jìn)一步構(gòu)造反演方程, 并利用l1范數(shù)最小化重構(gòu)稀疏信號(hào), 實(shí)現(xiàn)了對(duì)速度擾動(dòng)空間分布的定位.與線性最小二乘差分法相比, 稀疏重構(gòu)方法避免了復(fù)雜的參數(shù)確定操作, 克服了多擾動(dòng)點(diǎn)識(shí)別困難、定位不準(zhǔn)確等問題, 在保證反演精度的前提下能夠明顯減少計(jì)算時(shí)間.仿真結(jié)果表明, 稀疏重構(gòu)方法是一個(gè)性能優(yōu)越、穩(wěn)定可靠的優(yōu)秀算法, 有助于加快尾波干涉成像技術(shù)在地震火山預(yù)測(cè)、材料無損檢測(cè)等領(lǐng)域的實(shí)際應(yīng)用.

需要特別指出的是, 線性最小二乘法和稀疏重構(gòu)方法對(duì)于相對(duì)擾動(dòng)大小的定量分析效果均不佳,敏感核的質(zhì)量是主要影響因素.因此, 如何構(gòu)造能夠準(zhǔn)確描述尾波時(shí)空傳播特性的敏感核函數(shù)是尾波干涉成像的一個(gè)重要方面, 也是后續(xù)有待研究的問題.同時(shí), 本文將敏感核矩陣作為觀測(cè)矩陣, 其特性對(duì)稀疏重構(gòu)算法性能的影響有待后續(xù)著重研究.另外, 本文中的模擬實(shí)驗(yàn)在介質(zhì)尺寸與邊界條件等方面與數(shù)值仿真環(huán)境差距較大, 后續(xù)將設(shè)計(jì)實(shí)施更為合理的實(shí)驗(yàn)方案并進(jìn)一步研究多擾動(dòng)尾波干涉成像問題.