基于稀疏重構的尾波干涉成像方法*

張濤 侯宏? 鮑明

1) (西北工業大學航海學院, 海洋聲學信息感知工業和信息化部重點實驗室, 西安 710072)

2) (中國科學院噪聲與振動重點實驗室(聲學研究所), 北京 100190)

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

1 引 言

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

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

2 尾波干涉成像基本原理

2.1 尾波干涉基本原理

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

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

圖1 多散射介質中擾動前后波形的比較 (a) 直達波擾動前后的波形; (b)尾波擾動前后的波形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表示要分析的尾波部分的時間窗長度;t為時間窗中心位置;ts表示互相關中所用的走時差;R的值表示兩列波的相似程度.

當ts的某取值tsmax使得R(ts) 取最大值時, 該走時差tsmax即為所分析尾波部分擾動前后的時間延遲δt.那么根據Snieder[1]對尾波干涉的系統闡述, 可得兩列信號的相對波速變化

式中 δυ表示波速擾動, 是整體介質波速變化的加權平均;υ為擾動前的介質波速.

2.2 敏感核

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

式中的p(s1,s2,t) 表示波列從s1經過時間t到達s2的概率.在實際應用中, 此概率可以用波場強度等效, 波場強度是位置和時間的函數, 可以通過擴散方程的解來近似.在無限大二維介質中, 擴散方程的解為

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

2.3 反演成像模型

基于擴散近似敏感核, 文獻[10,11,13]提出了一個線性模型, 將尾波干涉法獲取的平均時延數據通過敏感核項K(s1,s2,x0,t) 與相應的局部變化聯系起來.

圖2 基于擴散近似的二維敏感核示例 (a) t= 1 s時的敏感核空間分布; (b) t= 5 s 時的敏感核空間分布; (c) t= 1 s 時的敏感核俯視圖; (d) t= 5 s時的敏感核俯視圖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.

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

式中,T∈RM,1表示通過不同的激勵接收對和不同的尾波觀察窗口獲取的M個時延數據;V∈RN,1表示空間各網格點對應的相對波速擾動值,N為網格總數;G∈RM,N表示T到V的映射,其矩陣元素Gij=?sKij, 其中Kij為敏感核矩陣K的元素,K∈RM,N的每行代表一個敏感核, 每列代表一個網格, 每個網格的面積為 ?s.

利用線性最小二乘差分法對V值進行反演, 即

式中,V0為先驗初始值, 一般為零矩陣;GT表示G的轉置;Cd為測得數據的對角協方差矩陣;Cm為介質空間元素平滑矩陣, 可用下式計算:

式中,σm為先驗標準差;λ0為網格間距;λ為相關長度.σm和λ一般通過 L 曲線法[25]選取.

在反演計算時, 還需要循環迭代Cm以尋找最優的V值, 迭代公式如下:

3 基于稀疏重構的尾波干涉成像方法

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

首先, 對信號進行稀疏表示如下:

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

其次, 構建一個與變換基底ψ不相關的觀測矩陣對a進行線性觀測

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

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

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

因此, (12)式的問題轉化為(13)式的凸優化問題,可以通過線性規劃理論求解, 本文采用應用廣泛的OMP法作為重構算法.

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

本文引入稀疏變換矩陣P對V進行離散稀疏變換.從而有

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

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

4 數值仿真與結果

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

本數值仿真運用二維聲波方程模擬波的傳播,

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

圖3 二維速度場模型Fig.3.2-D velocity field model.

圖4 激勵源及接收點布設Fig.4.Layout of the source and receivers.

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

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

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

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

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

圖7展示了實測數據的處理結果.該實驗是一個小型的模擬實驗, 利用激振器發出前述雷克子波, 作為對 1.5m×1.2m 的石膏板小樣品的激勵信號, 通過石膏板上布設的五個加速度傳感器采集微小擾動發生前后的波形信號.石膏板的波速與石膏板的含水率有關, 含水率增加, 板中超聲波傳播速度相應減小, 因此, 本實驗采用加濕的方式對樣品介質添加擾動.

將石膏板剖分為 1 0×10 的網格, 則待反演量的個數為100, 網格單元長寬分別為 0.15 m和0.12 m, 在石膏板的中心位置加水以增加中心點處的含水率作為局部擾動; 尾波分析時段為0.2—2.8 s,窗口長度均為 0.6 s, 重疊 0.4 s, 共 11 段尾波數據,則可獲得 5 ×11=55 個時延數據; 敏感核構造中的散射系數D=140m2/s ; 線性最小二乘法中的相對 長度λ=0.1m , 網格 間距λ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 實驗數據處理結果 (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.

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

前面討論了單點擾動算例, 并驗證了所提方法的優越性能.在實際中, 往往會出現多個擾動的情況, 因此算例3至5考慮了多個擾動區域的仿真,圖8—圖10分別給出了兩種算法反演結果的對比,圖中的藍色方框和黃色方框分別表示在該區域內施加+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次的結果, L曲線法獲 取 的σm分 別 為 0.01244, 0.00201和 0.00208,圖8(b), 圖9(b)和圖10(b)分別為算例 3, 4, 5 用本文方法迭代5次、5次和4次的結果.從圖9和圖10可以直觀地看出, 本文方法對于多擾動區域的反演質量明顯優于線性最小二乘法, 線性最小二乘法在算例4和算例5的結果中, 僅能定位出部分擾動, 且其位置與范圍也有一定程度的偏離真值,而本文方法對于設定的多擾動區域, 均能較為準確地反演出其位置與范圍; 而從算例3的反演圖像來看, 圖8(a)確定的擾動區域與真值相差甚遠, 甚至失去了參考價值, 相比之下, 圖8(b)的結果更為接近真值, 有利于正確分析引起波速變化的物理機制.

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

此外, 兩種方法在5個算例下反演成像時間對比如表1所列.

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

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

5 結 論

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

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