基于張量環分解的三維地震數據重建方法

張杏莉,劉作剛,張亞萍,趙衛東

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

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

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

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

1 方法原理

1.1 張量環分解

本研究采用Kolda等[21]和Long等[22]所提出的符號和定義。張量環分解是通過一系列低維核張量上的循環多線性乘積表示大維張量,這一系列低維核張量稱為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,張量環分解的定義為:

(1)

圖1 張量環分解

1.2 張量環低秩因子

將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通常在相應的數據中輸入觀測值更新,并且在每次迭代中通過TR因子[C]近似缺失的數據,即:

(6)

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

2.1 秩的選擇

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

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

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

2.2 TRLRF方法及復雜度分析

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

方法 1. 基于張量環分解的地震數據重建方法 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的階次比較低,但時間復雜度中含有數據維度平方相乘項,DDTF方法也包含維度平方相乘項,而TRLRF方法復雜度僅為維度平方之和,所以數據量較大時,TRLRF方法的時間復雜度比DDTF和OMPHR方法低。

3 實驗分析

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

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

3.1 實驗結果評價指標

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

(7)

(8)

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

3.2 仿真數據實驗及分析

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

圖3 原始仿真數據與地震道隨機缺失50%的仿真數據

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

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

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

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

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

圖6 仿真數據的不同采樣率與信噪比的關系

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

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

3.3 真實數據實驗及分析

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

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

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

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

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

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

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

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

圖11 真實數據的不同采樣率與信噪比的關系

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

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

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

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

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

圖12 地震道規則缺失50%的三維真實數據

圖13 地震道規則缺失50%的三維真實數據的各切面示例

圖14 地震道規則缺失50%的三維真實數據的側切面重建結果

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

表5 真實數據地震道規則缺失時重建結果對比

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

表6 真實數據地震道規則缺失時各切面重建結果對比

4 結論

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