李言言，李志明，王 琴

(中國地質(zhì)大學(xué) 數(shù)學(xué)與物理學(xué)院，湖北 武漢430074)

1 引 言

受采集環(huán)境和經(jīng)濟因素的影響，采集的地震數(shù)據(jù)常存在不規(guī)則分布，嚴重影響后續(xù)地震資料的處理和解釋。因此，對缺失地震數(shù)據(jù)進行插值在地震勘探中具有重要意義[1]。

目前，地震數(shù)據(jù)插值的方法主要有五類?；谙∈栊韵闰灥姆椒╗2]假設(shè)地震數(shù)據(jù)在變換域上具有稀疏性，將地震數(shù)據(jù)插值問題建模為對地震數(shù)據(jù)添加稀疏約束的數(shù)學(xué)反問題。常用的變換方法有Fourier變換[3]、Radon變換[4]、Curvelet變換[5]、Dreamlet變換[6]、Seislet變換[7]等?；陬A(yù)測濾波的方法[8]根據(jù)線性同相軸數(shù)據(jù)的可預(yù)測性對地震數(shù)據(jù)進行插值?；诓▌臃匠痰姆椒╗9，10]根據(jù)地下速度模型，利用偏移和反偏移算子對數(shù)據(jù)插值。近年來，基于深度學(xué)習(xí)的地震數(shù)據(jù)插值方法取得了一定的進展[11，12]。該方法通過端到端的學(xué)習(xí)缺失數(shù)據(jù)到完整數(shù)據(jù)的映射，并將訓(xùn)練好的網(wǎng)絡(luò)直接用于缺失數(shù)據(jù)插值。但是，該方法對數(shù)據(jù)的需求量較大，網(wǎng)絡(luò)泛化性不高。

基于秩減理論的插值方法近幾年引起了廣泛的關(guān)注[13-16]，這類方法假設(shè)完整的地震數(shù)據(jù)經(jīng)過預(yù)變換后具有低秩性，而地震數(shù)據(jù)的缺失會增加變換后數(shù)據(jù)的秩。因此，三維地震數(shù)據(jù)插值問題建模為矩陣或張量降秩問題。Oropeza等[15]理論上證明了線性同相軸的三維地震數(shù)據(jù)頻率切片所構(gòu)建的塊Hankel矩陣具有低秩性，提出將多道奇異譜分析方法用于地震數(shù)據(jù)的插值和去噪。Kreime等[17]基于高階奇異值分解(High-order Singular Value Decomposition, HOSVD)將其擴展到5D地震數(shù)據(jù)重建領(lǐng)域。由于構(gòu)建塊Hankel預(yù)變換的重建方法涉及較大矩陣的SVD(Singular Value Decomposition, SVD)分解，時間成本較高。Ma[18]提出紋理塊預(yù)變換方法，將三維數(shù)據(jù)體向量化得到二維矩陣，將地震插值問題轉(zhuǎn)化為矩陣降秩問題，采用APG (The Accelerated Proximal Gradient)和LMaFit (Low-rank Matrix Fitting)算法進行求解。Trickett等[19]提出Hankel張量預(yù)變換方法，基于頻率切片構(gòu)造四維Hankel張量，并采用低秩張量分解模型——PARAFAC(Parallel Factors)[20]求解，實現(xiàn)了較好的插值效果。張雪敏等[21]基于張量Tucker秩分解方法對張量降秩，利用并行矩陣分解的低秩張量補全 (Low-rank Tensor Completion by Parallel Matrix Factorization, Tmac)算法，無SVD分解步驟，提高了重建效率。但是，上述兩種方法都需要預(yù)先估計張量的秩。Liu等[22]將插值問題建模為張量秩最小化模型，基于張量Tucker分解，將張量秩凸松弛為張量核范數(shù)求解，計算簡單，無需預(yù)先估計張量秩。但是該方法忽略張量每一個方向模展開矩陣的奇異值的先驗信息[23]，僅最小化各模展開矩陣的奇異值之和，得到原始秩問題的次優(yōu)解。

研究表明，相比于核范數(shù)，非凸函數(shù)更近似矩陣秩，而且已經(jīng)在圖像修復(fù)、信號處理和地震數(shù)據(jù)重建領(lǐng)域取得了較好的應(yīng)用效果[24]。Zhang等[25，26]提出非凸log-sum函數(shù)和截斷核范數(shù)正則化方法用于二維地震數(shù)據(jù)插值。相比于凸松弛方法，非凸模型明顯地提高了地震數(shù)據(jù)插值精度。本文提出張量秩的非凸替代，利用logε函數(shù)代替張量秩，并將其轉(zhuǎn)化為張量的加權(quán)核范數(shù)模型，根據(jù)加權(quán)閾值算子得到其最優(yōu)閉式解。為驗證本文方法，與基于塊Hankel預(yù)變換的經(jīng)典MSSA(Multichannel Singular Spectrum Analysis)方法和基于Hankel張量預(yù)變換的LRTC(Low-rank Tensor Completion)方法作對比。實驗結(jié)果說明，本文方法實現(xiàn)了更高的插值精度。

2 方 法

2.1 張量的基本知識及符號表示

張量有多種秩的定義[20]，其中，張量Tucker秩以矩陣秩為基礎(chǔ)，易于求解，得到了廣泛的應(yīng)用。本文基于張量Tucker秩建立地震數(shù)據(jù)插值的數(shù)學(xué)模型。

首先介紹張量的相關(guān)符號表示。記r階張量X∈RI1×I2×…×Ir每一個元素為Xj1j2...jr=X(j1,j2,...,jr),其中ji={1,...,Ii}。固定任意r-1個維度，可得張量的纖維向量，如X (j1,j2,...,jr-1,:)∈RIr×1。X的Tucker秩的定義基于X的模展開[20]：

(1)

其中，rank(·)是矩陣的秩函數(shù)；X(i)∈RIi×(I1×…Ii-1×Ii+1×…Ir),其中(i=1,2,...,r)是X的模-i展開矩陣；X(i)的秩是X的Tucker秩的第i個分量。X與X(i)的元素對應(yīng)關(guān)系為：

Xj1j2...jr=X(i)(ji,m)

(2)

其中,

記張量的模展開算子為：

X(i)=unfoldi(X)

(3)

記反變換為X=foldi(X(i))。

2.2 地震數(shù)據(jù)的Hankel張量預(yù)變換

完整地震數(shù)據(jù)構(gòu)建的Hankel張量具有低秩性[19]，本文基于Hankel張量預(yù)變換對缺失地震數(shù)據(jù)進行插值。

考慮三維地震記錄Dt=D(t,x,y)∈RNt×Nx×Ny。其中，Nt表示時間采樣點數(shù)目；Nx,Ny表示沿空間方向x,y的地震道數(shù)目。沿Dt的時間方向作一維Fourier變換，得到頻域數(shù)據(jù)Dω=D(ω,x,y)∈CNω×Nx×Ny，固定頻率ω=ωj，得到頻率切片Dωj=D(ωj,x,y)∈CNx×Ny。

根據(jù)Dωj構(gòu)建Hankel張量：

Yωj=H(Dωj)

(4)

圖1 Hankel張量變換示意圖Fig.1 Schematic diagram of Hankel tensor transformation

2.3 地震數(shù)據(jù)插值的數(shù)學(xué)模型

基于Hankel張量的低秩性，可以通過對Hankel張量Yωj降秩實現(xiàn)地震數(shù)據(jù)插值。建立式(5)所示的地震數(shù)據(jù)3D插值的數(shù)學(xué)模型：

(5)

其中，ranktc(·)是張量的Tucker秩;Xωj是待重建的四維張量;Yωj是Hankel張量數(shù)據(jù);Ω是Yωj的非零元素對應(yīng)的的索引指標(biāo)集;PΩ是采樣算子。為方便表示，下文省略張量的下標(biāo)ωj。

根據(jù)張量Tucker秩的定義，式(5)是多目標(biāo)優(yōu)化問題，一般地，將張量的Tucker秩轉(zhuǎn)化為模-i展開矩陣秩的加權(quán)和，見式(6)：

(6)

(7)

事實上，奇異值包含數(shù)據(jù)的先驗信息，較大奇異值包含數(shù)據(jù)的主要信息，較小奇異值表示數(shù)據(jù)的噪聲。然而，核范數(shù)最小化的同時最小化了所有奇異值之和，不能保護數(shù)據(jù)的有效信息，使得基于該方法的重建精度得不到保證。本文提出非凸張量秩最小化(Nonconvex Tensor Rank Minimization, NTRM)模型，用非凸logε模型代替核范數(shù)，并采用迭代軟閾值算子自適應(yīng)對奇異值加權(quán)，改善重建質(zhì)量。所提出的NTRM模型見式(8)：

(8)

圖2展示了秩為1時，核范數(shù)、rank函數(shù)和非凸logε函數(shù)(基于不同的ε值)隨奇異值變化的情況。由圖2可知，相比于核范數(shù)，非凸logε函數(shù)隨奇異值變化較小，擬合秩呈現(xiàn)較好的效果。

圖2 核范數(shù)、真實秩和logε函數(shù)隨奇異值變化曲線Fig.2 Nuclear norm, true rank and logεfunction with singular value curve

2.4 NTRM模型求解

對于式(8)的非凸模型，首先，引入中間變量M(i)(其中i=1,2,3,4)，將式(8)轉(zhuǎn)化為式(9)求解：

X(i)=M(i)

(9)

構(gòu)建式(9)的增廣拉格朗日函數(shù)，將式(9)轉(zhuǎn)化為：

(10)

其中，Λ(i)是拉格朗日乘子。本文利用ADMM (Alternating Direction Method of Multipliers)框架迭代求解X,M(i),Λ(i),λ，更新式如下(n表示迭代次數(shù))：

(11)

(12)

(13)

λ(n+1)=tλ(n)(t=1.05)

(14)

(15)

根據(jù)式(12)，易求得X(n+1)的閉式解：

(16)

表1 算法1:基于非凸張量秩最小化(NTRM)方法的地震數(shù)據(jù)插值Table 1 Algorithm 1: seismic data interpolation based on nonconvex tensor rank minimization method

3 數(shù)據(jù)實驗

本文通過仿真和真實疊前數(shù)據(jù)驗證NTRM模型的有效性，并與傳統(tǒng)的MSSA方法和LRTC方法作對比。使用信噪比(Signal-to-Noise, SNR)作為評價指標(biāo)：

(17)

3.1 仿真數(shù)據(jù)實驗

仿真地震數(shù)據(jù)包含3個彎曲同相軸，101個時間采樣點，x,y方向各包含32個地震道。圖3(a)表示原始完整地震數(shù)據(jù);圖3(e)表示隨機缺失70 %的地震數(shù)據(jù),該缺失地震數(shù)據(jù)的SNR值為1.55 dB。在圖3(b)～ 圖3(d)中展示的是基于MSSA、LRTC、NTRM方法的重建結(jié)果圖。其中，對于經(jīng)典的MSSA方法選取秩為10，此時SNR值最高，為29.99 dB。LRTC和NTRM方法均基于Hankel張量預(yù)變換，設(shè)置各模-i展開矩陣的權(quán)值向量α=[0.25,0.25,0.25,0.25]，設(shè)置NTRM模型中的logε函數(shù)的超參數(shù)ε=0.001，重建SNR值分別是32.43 dB和55.46 dB。為更直觀地說明NTRM模型的優(yōu)勢，在圖3(f)～ 圖3(h)中展示各個方法對應(yīng)的重建殘差圖像。由重建結(jié)果圖和殘差圖像可知，NTRM方法的重建結(jié)果最近似于原始地震圖像，具有較高的重建質(zhì)量。

表1展示了不同缺失率下，三種方法的重建結(jié)果的SNR值對比情況。由表2可知，隨著缺失率的增加，數(shù)據(jù)重建的SNR減小。當(dāng)缺失率為80 %時，LRTC方法和MSSA方法的重建精度較低，但是，NTRM方法依然達到較高的重建精度。

為了驗證本文方法的抗噪性能，在POCS(Projection onto Convex Sets)框架下同時完成插值和去噪[15]，具體迭代式如下：

(18)

其中，an(n是迭代次數(shù))是與迭代次數(shù)相關(guān)的標(biāo)量，滿足a1=1,anmax=0(nmax是總迭代次數(shù))。Ωc是Ω的補集。本文通過仿真數(shù)據(jù)測試去噪效果。首先，對仿真數(shù)據(jù)添加均值為0，標(biāo)準(zhǔn)差為0.01的高斯噪聲(圖4a)。隨機缺失50 %地震道的數(shù)據(jù)如圖4(b)所示。缺失數(shù)據(jù)的信噪比為2.76 dB。在圖4(c)和圖4(d)中展示MSSA方法(基于POCS框架同時插值和去噪)和NTRM的重建結(jié)果。重建信噪比分別是23.96 dB和26.16 dB。由重建圖像和重建信噪比可知，本文所提方法也具有較好的抗噪性能。

圖4 基于MSSA方法和NTRM方法的重建結(jié)果對比Fig.4 Comparison of reconstruction results based on MSSA method and NTRM method

3.2 真實數(shù)據(jù)實驗

如圖5(a)所示的完整地震數(shù)據(jù)包含200個時間采樣點，采樣間隔是4 ms，x,y方向分別包含31，10地震道。圖5(e)表示隨機缺失60 %的地震數(shù)據(jù)。此缺失數(shù)據(jù)的SNR值為2.46 dB?；贛SSA、LRTC 和NTRM方法的重建結(jié)果圖分別在圖5(b)～ 圖5(d)中展示，重建信噪比分別為16.45 dB、12.71 dB、19.14 dB。由圖5(b)～ 圖5(d)的重建結(jié)果和重建信噪比可知，相比于展示的另外兩種方法，基于NTRM模型的重建精度最高。圖5(f)～ 圖5(h)中展示的是MSSA、LRTC和NTRM方法重建結(jié)果對應(yīng)的殘差圖像。由圖5(f)～圖5(h)可知，NTRM方法對應(yīng)的殘差具有較小的幅值，重建數(shù)據(jù)最接近原始數(shù)據(jù)。圖6展示了缺失率分別為30 %、40 %、50 %、60 %、70 %時，MSSA、LRTC、NTRM方法的重建SNR值對比情況，可以看出，即使在高缺失率時，NTRM模型也具有較好的重建效果。

圖5 基于不同方法的地震數(shù)據(jù)重建結(jié)果對比Fig.5 Comparison of seismic data reconstruction results based on different methods

圖6 不同缺失率時，基于MSSA、LRTC和NTRM方法的三維缺失地震數(shù)據(jù)重建SNR值對比Fig.6 Comparison of SNR values reconstructed from 3D seismic data with different miss rates based on MSSA, LRTC and NTRM methods

4 結(jié) 論

經(jīng)典的MSSA方法通過直接截斷塊Hankel矩陣奇異值用于地震數(shù)據(jù)插值，得到的解往往不是最優(yōu)的。地震數(shù)據(jù)的Hankel張量預(yù)變換方法重排塊Hankel為四維Hankel張量，具有較好的空間幾何結(jié)構(gòu)，張量低秩性得到保證?；贖ankel張量秩最小化的地震插值方法得到廣泛應(yīng)用。傳統(tǒng)的LRTC方法解決了張量秩難以求解的問題，但是，該方法用核范數(shù)近似秩，誤差較大。本文提出的NTRM模型，將二維非凸模型擴展到四維張量領(lǐng)域用于地震數(shù)據(jù)插值，提高了模型建立的準(zhǔn)確性，進而提高了插值質(zhì)量。對于非凸模型，本文利用加權(quán)軟閾值算子迭代求解。仿真和真實三維地震數(shù)據(jù)插值實驗表明，相比于經(jīng)典的MSSA和LRTC方法，NTRM方法的重建精度最高。但是本文方法僅用于隨機缺失地震數(shù)據(jù)插值，對規(guī)則缺失數(shù)據(jù)插值不適用。今后，筆者將把Cadzow重建理論與本文方法相結(jié)合用于規(guī)則缺失道地震數(shù)據(jù)插值作為下一步的研究計劃。