張杏莉,劉作剛,張亞萍,趙衛(wèi)東

(山東科技大學 計算機科學與工程學院,山東 青島 266590)

地震數(shù)據(jù)采集受采集環(huán)境和地理條件的限制,容易造成數(shù)據(jù)不完整,低質(zhì)量的地震數(shù)據(jù)對后續(xù)地震數(shù)據(jù)處理和解釋工作產(chǎn)生較大影響,因此地震數(shù)據(jù)重建已成為地震數(shù)據(jù)處理的重點研究內(nèi)容之一。數(shù)據(jù)重建是在數(shù)據(jù)缺失的情況下,利用少量的數(shù)據(jù)預測、補全缺失的數(shù)據(jù)。目前常用的數(shù)據(jù)重建方法有濾波器、函數(shù)變換、壓縮感知和降秩等。

基于濾波器的方法通過缺失數(shù)據(jù)與某種插值濾波器的卷積來獲得重建數(shù)據(jù)。Spitz[1]提出f-x域插值重建方法,在f-x域中使用低頻數(shù)據(jù)重建高頻數(shù)據(jù);Naghizadeh等[2]提出線性預測濾波器重建方法,在規(guī)則的網(wǎng)格下重建不規(guī)則的數(shù)據(jù)。但基于濾波器的方法只能處理規(guī)則的數(shù)據(jù)。基于函數(shù)變換的方法先對地震數(shù)據(jù)進行正變換編碼,然后在變換域內(nèi)插值,再進行反變換解碼,進而對數(shù)據(jù)重建。Jahanjoo等[3]提出抗泄漏傅里葉變換方法,通過正則化處理不規(guī)則的數(shù)據(jù),但是當輸入的數(shù)據(jù)比較混疊時數(shù)據(jù)重建效果不理想,并且傅里葉變換是一種全局變換,不能準確地表示局部特征;唐歡歡等[4]提出3D高階拋物Radon變換方法,解決了3D不規(guī)則數(shù)據(jù)的重建問題,但Radon變換處理復雜曲線事件的能力較差;Zhang等[5]通過多尺度、多方向的Curvelet變換對較大規(guī)模的數(shù)據(jù)進行處理;文獻[6]采用一種非等間距的離散Curvelet變換方法對非均勻采樣的數(shù)據(jù)進行處理;張凱等[7]將Shearlet變換和離散余弦變換字典相結合進行數(shù)據(jù)重建。但是函數(shù)變換方法的基函數(shù)是固定的,不能根據(jù)地震數(shù)據(jù)的特征自適應地選擇基函數(shù)[8]。壓縮感知理論是一種新的信號處理方法,可以由遠少于采樣定理要求的采樣點重建恢復,目前已經(jīng)被應用到地震數(shù)據(jù)重建中。Zhang等[9]提出基于隨機樣本的壓縮感知方法,通過將樣本劃分為若干子集以獲得重建結果;孫苗苗等[10]將形態(tài)成分分析與壓縮感知結合,消除了空間假頻;段中鈺等[11]將壓縮感知和平方正則交替方向乘子法結合,消除了過擬合現(xiàn)象。但壓縮感知方法難以滿足地震道大范圍缺失的情況。降秩是地震數(shù)據(jù)重建的另一種方法,如低秩矩陣補全方法[12]、基于奇異值分解(singular value decomposition, SVD)的低秩約束和曲線域稀疏約束的重建方法[13-14],而在低秩矩陣上進行核范數(shù)約束的方法變得越來越廣泛。Wang等[15]提出加權核范數(shù)最小化方法,通過迭代重新加權使加權核函數(shù)模型接近于“計數(shù)秩函數(shù)”并且可以提升低秩矩陣,但核范數(shù)最小化是指所有奇異值之和的最小化,不等于核函數(shù)的最小值;Zhang等[16]對核范數(shù)正則化,從而對更少的奇異值最小化,可以更加逼近秩函數(shù),但是奇異值分解會消耗大量的計算資源;多通道奇異譜分析方法是一種數(shù)據(jù)驅(qū)動方法,該方法基于Hankel矩陣的截斷奇異值分解,通過頻率切片進行數(shù)據(jù)重建[17]。但是傳統(tǒng)的低秩或降秩矩陣補全方法在三維數(shù)據(jù)的應用中仍存在一些重建誤差。

張量環(huán)(tensor ring, TR)分解是一種新的數(shù)據(jù)補全方法,對數(shù)據(jù)特征沒有要求,可以處理復雜的事件,并且可以有效利用三維數(shù)據(jù)的結構信息,被廣泛應用于圖像的缺失數(shù)據(jù)重建,如子空間非局部張量環(huán)分解的高光譜圖像重建[18]、全變分正則化張量環(huán)分解的遙感圖像重建[19]、張量環(huán)低秩因子(tensor ring low-rank factors, TRLRF)缺失圖像重建[20]等。本研究將TRLRF方法首次應用于地震數(shù)據(jù)重建,先通過張量環(huán)分解將三維地震數(shù)據(jù)轉(zhuǎn)換為小的三維數(shù)據(jù)的乘積,再利用交替方向乘子法(alternating direction method of multipliers, ADMM)和增廣拉格朗日函數(shù)對小的三維數(shù)據(jù)求解,最后通過循環(huán)多線性乘積將小的三維數(shù)據(jù)恢復為大的三維數(shù)據(jù),實現(xiàn)對三維地震數(shù)據(jù)的重建。仿真數(shù)據(jù)和真實數(shù)據(jù)實驗結果均表明,TRLRF方法在重建效率和重建質(zhì)量上都有較大優(yōu)勢,重建效果良好。

1 方法原理

1.1 張量環(huán)分解

本研究采用Kolda等[21]和Long等[22]所提出的符號和定義。張量環(huán)分解是通過一系列低維核張量上的循環(huán)多線性乘積表示大維張量,這一系列低維核張量稱為TR因子。對于n=1,2,…,N,TR因子可以用C(n)∈RRn×In×Rn+1表示,其中Rn(1≤n≤N)表示TR的秩,R1=Rn+1,In是原始張量的維度。對于一個

N階張量A∈RI1×I2×…×IN,張量環(huán)分解的定義為:

(1)

圖1 張量環(huán)分解

1.2 張量環(huán)低秩因子

將ADMM和增廣的拉格朗日方程相結合,TRLRF方法的表達式為[20]:

(2)

C(n)、W(n,i)和U(n,i)都是相對獨立的,式(2)可以通過以下過程更新。

1) 更新C(n)。

(3)

2) 更新W(n,i)。

(4)

3) 更新U(n,i)。

U(n,i)=U(n,i)+η(W(n,i)-C(n))。

(5)

4) 更新A。A通常在相應的數(shù)據(jù)中輸入觀測值更新,并且在每次迭代中通過TR因子[C]近似缺失的數(shù)據(jù),即:

(6)

2 TRLRF方法的影響因素及復雜度分析

2.1 秩的選擇

在張量環(huán)分解中一個關鍵問題是如何選取最優(yōu)秩。一個地震數(shù)據(jù)體包含的地震信息都是相互關聯(lián)的,秩選取過小會導致TR因子中的數(shù)據(jù)信息過少,秩選取過大會導致TR因子中存在過多的干擾數(shù)據(jù)。因此,秩選取過大或過小均會使TR因子不能準確地近似缺失的數(shù)據(jù),導致重建效果差。同時,秩選取過大會增加方法運行時間,影響重建效率。本研究通過對秩取不同值時的信噪比(signal-to-noise ratio, SNR)、均方根誤差(root mean square error, RMSE)以及運行時間等評價指標進行統(tǒng)計分析確定最優(yōu)秩的值。

分別使用兩個包含線性事件(數(shù)據(jù)1、數(shù)據(jù)2)、兩個包含曲線事件(數(shù)據(jù)3、數(shù)據(jù)4)的仿真數(shù)據(jù)對秩的選擇進行實驗,并計算重建后數(shù)據(jù)的SNR、RMSE和運行時間,實驗結果如圖2所示。由圖2(a)可知,隨著秩的增加,重建后的仿真數(shù)據(jù)的SNR值逐步升高,當秩為10時SNR值取得最大值,當秩大于10時SNR值趨于穩(wěn)定或稍有下降趨勢;由圖2(b)可知,隨著秩的增加,RMSE值逐漸降低,當秩為10時RMSE取得最小值,當秩大于10時RMSE趨于平穩(wěn)或逐漸增加?？梢?秩為10時4個仿真數(shù)據(jù)的重建效果最好。由圖2(c)可知,隨著秩的增加,TRLRF方法所需運行時間快速增加。因此,綜合考慮重建后數(shù)據(jù)的SNR、RMSE以及重建所需的運行時間,認為秩的值取10時可以獲得最好的重建結果并且保持良好的重建效率。因此,TRLRF方法中秩的值選擇為10。

圖2 不同秩時的SNR、RMSE和運行時間

2.2 TRLRF方法及復雜度分析

表1 不同重建方法的時間復雜度

方法 1. 基于張量環(huán)分解的地震數(shù)據(jù)重建方法 1) Input:missing seismic data PΩ(T), TR-rank{Rn } N n= 1 = 10 2) Initialization:For n = 1,…,N,i = 1,2,3, random sample C (n) by N ～ (0,1),ξ = 10,η = 0. 9,tol = 10 -8 ,kmax = 500,W (n,i) = 0,U (n,i) = 0. 3) For k = 1 to kmax do 4) Apre = A . 5) Calculate (C (≠n),T (2) C (≠n) (2) ) -1 ,A < n > C (≠n) (2) . 6) Update {C (n) } N n = 1 by (3). 7) Calculate S 1 η C (n) (i) - 1 η U (n,i) (i) ) . 8) Update {W (n,i) } N,3 n = 1,i = 1 by (4). 9) Update {U (n,i) } N,3 n = 1,i = 1 by (5). 10) Update A by (6). 11) If mse(A - Apre) < tol ,break. 12) End for 13) Output:reconstructed seismic data

由表1可見,盡管OMPHR方法中K的階次比較低,但時間復雜度中含有數(shù)據(jù)維度平方相乘項,DDTF方法也包含維度平方相乘項,而TRLRF方法復雜度僅為維度平方之和,所以數(shù)據(jù)量較大時,TRLRF方法的時間復雜度比DDTF和OMPHR方法低。

3 實驗分析

為驗證TRLRF方法的有效性,分別對仿真數(shù)據(jù)的地震道隨機缺失和實際三維地震數(shù)據(jù)(來源https:∥wiki.seg.org/wiki/Kerry-3D)的地震道隨機缺失、規(guī)則缺失進行實驗。

在地震數(shù)據(jù)重建方法中,把數(shù)據(jù)劃分為字典小塊的字典學習和對數(shù)據(jù)的每個二維頻率切片的矩陣補全的重建方法都是典型方法。近年來,有學者將深度學習應用于地震數(shù)據(jù)重建,但是需要大量樣本數(shù)據(jù)進行訓練和測試。故本研究選用基于字典學習的DDTF方法和基于矩陣補全的OMPHR方法作為本研究的對比方法。

3.1 實驗結果評價指標

本研究采用地震數(shù)據(jù)重建后的SNR、RMSE和運行時間3個指標驗證數(shù)據(jù)重建的效果和效率。SNR(SNR)和RMSE(RMSE)兩種指標分別定義為:

(7)

(8)

式中:f是完整原始地震數(shù)據(jù),f0是重建地震數(shù)據(jù),mean(·)表示元素的平均值。信噪比越高表明重建數(shù)據(jù)越接近原始數(shù)據(jù),均方根誤差是重建后的地震數(shù)據(jù)偏離原始數(shù)據(jù)平均值的度量,均方根誤差越小表示重建的數(shù)據(jù)越接近原始數(shù)據(jù)。

3.2 仿真數(shù)據(jù)實驗及分析

為驗證TRLRF方法的有效性和優(yōu)越性,使用Marmousi2構造仿真地震數(shù)據(jù),圖3(a)為原始仿真數(shù)據(jù)的正面切片,圖3(b)為隨機缺失50%的仿真數(shù)據(jù)。從圖3(b)可以看出,在缺失50%的數(shù)據(jù)后,該地震數(shù)據(jù)損壞嚴重。分別使用TRLRF、DDTF和OMPHR方法對圖3(b)所示地震數(shù)據(jù)進行重建,重建后的地震數(shù)據(jù)如圖4所示。

圖3 原始仿真數(shù)據(jù)與地震道隨機缺失50%的仿真數(shù)據(jù)

圖4 地震道隨機缺失50%的仿真數(shù)據(jù)的正切面重建結果

由圖4可知,原始仿真數(shù)據(jù)在采樣點600之前地震波的幅值較大,特征比較明顯,而在采樣點600之后地震波的幅值較小,波形趨于平緩。DDTF方法在重建過程中通過小塊數(shù)據(jù)形成訓練集恢復缺失的數(shù)據(jù),在幅值較小、特征不明顯的情況下恢復數(shù)據(jù)能力較差,故DDTF方法在采樣點600之后的重建效果較差。OMPHR方法通過二維切面的頻率切面和對每個切面進行反平均生成每一列以恢復缺失的數(shù)據(jù),在進行反平均生成每一列時容易造成數(shù)據(jù)的混亂,因此,OMPHR方法在重建后產(chǎn)生反向波的混亂現(xiàn)象。而TRLRF方法通過張量環(huán)分解有效利用了三維數(shù)據(jù)的結構信息,能保證劃分的塊包含更多的特征,使用適合解決結構復雜問題的交替方向乘子法求解劃分的塊,從而獲得更好的重建結果。

從微觀角度進一步對比3種方法的重建效果,對第65條缺失的單個信道的原始數(shù)據(jù)與重建結果進行對比,如圖5所示。由圖5可看出,TRLRF方法較完整地重建了缺失的數(shù)據(jù),而DDTF方法對采樣點520后幅值較小的波形重建效果差,OMPHR方法在采樣點500后有明顯振蕩,重建質(zhì)量較差。

圖5 不同方法重建后的第65條信道

為了驗證不同缺失率下3種方法的重建效果,對3種方法的SNR進行分析,結果如圖6所示。由圖6可以看出,TRLRF方法在不同缺失率下的重建質(zhì)量均優(yōu)于DDTF和OMPHR方法。

圖6 仿真數(shù)據(jù)的不同采樣率與信噪比的關系

分別計算TRLRF、DDTF和OMPHR方法在地震道隨機缺失50%重建后的SNR、RMSE和運行時間,計算結果見表2。分析表2可知,TRLRF方法在缺失率相等的條件下取得最高的SNR和最低的RMSE。由此可以看出,TRLRF方法在仿真數(shù)據(jù)的重建質(zhì)量上比DDTF和OMPHR方法都要優(yōu)越,并且TRLRF方法的運行時間最小,具有較高的重建效率。

表2 仿真數(shù)據(jù)地震道隨機缺失時重建結果對比

3.3 真實數(shù)據(jù)實驗及分析

使用真實地震數(shù)據(jù)(來源https:∥wiki.seg.org/wiki/Kerry-3D)進行實驗,如圖7、圖8所示。

圖7 原始三維真實數(shù)據(jù)與地震道隨機缺失50%的三維真實數(shù)據(jù)

圖8 原始三維真實數(shù)據(jù)與地震道隨機缺失50%的三維數(shù)據(jù)的各切面示例

分別使用TRLRF、DDTF和OMPHR方法對圖7(b)缺失的數(shù)據(jù)進行重建,其正切面的重建結果如圖9所示。由圖9可以看出,DDTF方法重建的數(shù)據(jù)存在大面積誤差,重建效果差;OMPHR方法重建的數(shù)據(jù)在某些地震道上存在誤差,重建效果較差;而TRLRF方法在整體重建質(zhì)量上優(yōu)于DDTF和OMPHR方法。這是由于DDTF方法是把數(shù)據(jù)劃分為小塊數(shù)據(jù)形成訓練集,通過訓練集對缺失的數(shù)據(jù)進行估計,訓練過程中對特征不明顯的數(shù)據(jù)不敏感,不能較好地補全缺失的數(shù)據(jù)。OMPHR方法是把三維數(shù)據(jù)的每個二維切面轉(zhuǎn)換為頻率切面,對每個切面采用OR1MP方法將矩陣表示為秩為1的基矩陣的線性組合,但是對于秩比較大的缺失數(shù)據(jù)來說,秩為1的基矩陣的線性組合不能很好地表示原始數(shù)據(jù),從而不能很好地恢復缺失的數(shù)據(jù)。而TRLRF方法通過張量環(huán)分解有效利用了三維數(shù)據(jù)的結構信息,保證劃分的塊包含更多地特征,使用適合解決結構復雜問題的交替方向乘子法求解劃分的塊,從而獲得更好的重建結果。

圖9 地震道隨機缺失50%的三維真實數(shù)據(jù)的正切面重建結果

從微觀角度進一步對比3種方法的重建效果,對正切面的第33條缺失的單個信道的原始數(shù)據(jù)與重建結果進行對比,如圖10所示。從圖10可以看出,TRLRF方法較完整的重建了缺失的數(shù)據(jù),而DDTF和OMPHR方法存在振蕩,重建質(zhì)量較差。

圖10 不同方法重建后正切面的第33條信道

為了說明3種方法在不同缺失率下的良好重建性能,對3種方法在不同缺失率下重建后的SNR進行分析,如圖11所示。由圖11可以看出,TRLRF方法在不同缺失率下的重建質(zhì)量均優(yōu)于DDTF和OMPHR方法。

圖11 真實數(shù)據(jù)的不同采樣率與信噪比的關系

分別計算TRLRF、DDTF和OMPHR方法在地震道隨機缺失50%重建后的SNR、RMSE和運行時間,計算結果見表3。分析表3可知,TRLRF方法的SNR最高,而OMPHR方法重建的SNR僅為TRLRF方法的1/2;TRLRF方法的RMSE值最小,可見,TRLRF方法重建后的三維地震數(shù)據(jù)最接近原始三維地震數(shù)據(jù);同時,TRLRF方法的運行時間相對于DDTF和OMPHR方法最短。因此,TRLRF方法具有良好的重建效率。

表3 真實數(shù)據(jù)地震道隨機缺失時重建結果對比

表4分別給出了真實數(shù)據(jù)地震道隨機缺失時上切面、側(cè)切面和正切面的SNR和RMSE結果。從表4可知,TRLRF方法在各個切面的重建效果優(yōu)于DDTF和OMPHR方法。

表4 真實數(shù)據(jù)地震道隨機缺失時各切面重建結果對比

對地震數(shù)據(jù)的規(guī)則缺失進行實驗,進一步說明TRLRF方法能夠很好地重建缺失數(shù)據(jù),如圖12、圖13所示。分別使用TRLRF、DDTF和OMPHR方法對圖12缺失的數(shù)據(jù)進行重建,其側(cè)切面的重建結果如圖14所示。由圖14可以看出,DDTF方法重建的數(shù)據(jù)存在大面積誤差,重建效果差;OMPHR方法重建的數(shù)據(jù)在某些地震道上存在誤差,重建效果較差;而TRLRF方法在整體重建質(zhì)量上優(yōu)于DDTF和OMPHR方法。這是由于DDTF方法是把數(shù)據(jù)劃分為小塊數(shù)據(jù)形成訓練集,通過訓練集對缺失的數(shù)據(jù)進行估計,在訓練過程中對特征不明顯的數(shù)據(jù)不敏感,不能較好地補全缺失的數(shù)據(jù)。OMPHR方法是把三維數(shù)據(jù)的每個二維切面轉(zhuǎn)換為頻率切面,難以充分利用三維數(shù)據(jù)的結構信息,不能很好地恢復缺失的數(shù)據(jù)。而TRLRF方法通過張量環(huán)分解有效利用了三維數(shù)據(jù)的結構信息,可保證劃分的塊包含更多特征,使用適合解決結構復雜問題的交替方向乘子法求解劃分的塊,從而獲得更好的重建結果。

圖12 地震道規(guī)則缺失50%的三維真實數(shù)據(jù)

圖13 地震道規(guī)則缺失50%的三維真實數(shù)據(jù)的各切面示例

圖14 地震道規(guī)則缺失50%的三維真實數(shù)據(jù)的側(cè)切面重建結果

分別計算TRLRF、DDTF和OMPHR方法在地震道規(guī)則缺失50%重建后的SNR、RMSE和運行時間,結果見表5。分析表5可知,TRLRF方法的SNR最高,是DDTF和OMPHR方法重建的SNR值的兩倍多;TRLRF方法的RMSE值最小,重建后的三維地震數(shù)據(jù)更接近原始三維真實地震數(shù)據(jù)。同時,TRLRF方法的運行時間均小于DDTF和OMPHR方法?？梢?TRLRF方法具有良好的重建效率。

表5 真實數(shù)據(jù)地震道規(guī)則缺失時重建結果對比

表6分別給出了真實數(shù)據(jù)地震道規(guī)則缺失時上切面、側(cè)切面和正切面的SNR和RMSE結果,可以看出TRLRF方法在各個切面上的重建效果均優(yōu)于DDTF和OMPHR方法。

表6 真實數(shù)據(jù)地震道規(guī)則缺失時各切面重建結果對比

4 結論

本研究針對三維地震數(shù)據(jù)缺失重建問題,提出一種基于隱空間上張量環(huán)低秩因子的地震數(shù)據(jù)重建方法。張量環(huán)分解方法對數(shù)據(jù)特征要求低,能夠處理幅值較低、特征不明顯的數(shù)據(jù),并且能夠充分利用三維數(shù)據(jù)的結構信息。仿真數(shù)據(jù)和真實地震數(shù)據(jù)的實驗結果表明,在SNR、RMSE和運行時間等量化指標下,與DDTF、OMPHR方法相比,TRLRF方法重建效果較好,同時針對較大的三維地震數(shù)據(jù),張量環(huán)分解可有效提高計算效率。本研究依據(jù)仿真三維地震數(shù)據(jù)在不同秩時的重建效果確定張量環(huán)分解中最優(yōu)秩,該方法具有一定的參考意義。復雜三維地震數(shù)據(jù)在進行張量環(huán)分解時的最優(yōu)秩的自適應選取是下一步的研究方向和研究重點。