非凸Lp范數(shù)三維地震數(shù)據(jù)重建

劉 群 付麗華 張婉娟

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

0 引言

地震數(shù)據(jù)為探測地質(zhì)構造和地層分布提供了有效的地球物理信息。地震數(shù)據(jù)的質(zhì)量直接影響地震處理和解釋結果[1-4]。地震數(shù)據(jù)在空間上可能出現(xiàn)不規(guī)則采樣[5-6]，造成有效信號的缺失，從而增加地震數(shù)據(jù)處理難度[7-8]。因此，對地震數(shù)據(jù)進行有效重建，使后續(xù)數(shù)據(jù)處理能夠更好地刻畫復雜地質(zhì)結構非常必要[9]。隨著地震勘探研究的不斷深入,更多的數(shù)學方法被用于地震數(shù)據(jù)重建[10-12]。其中一類是基于濾波的重建方法，該類方法通過褶積插值濾波器實現(xiàn)[13-16]。另一類重要方法是基于變換域的地震數(shù)據(jù)重建，主要有Radon域變換[17-19]、Fourier域變換[20-22]、Curvelet域變換[23-25]和Shearlet變換[26]等，即基于數(shù)學變換理論，把數(shù)據(jù)轉(zhuǎn)換到不同的域進行重建?；诓▌臃匠痰牡卣饠?shù)據(jù)重建方法也得到廣泛研究，如炮域延拓法[27]、炮檢距域延拓法[28]、DMO法[29]和拓展成像條件法[30]，該類方法計算量大，運算耗時，且當?shù)叵滦畔⑽粗蚓容^低時會影響重建效果，所以實際應用較少。

近年來，基于降秩重建的方法引起廣泛關注。Trickett等[31]和Naghizadeh等[32]提出基于Cadzow濾波的三維地震數(shù)據(jù)隨機噪聲衰減方法，建立了Hankel矩陣并對其進行奇異值分解，通過降秩壓制線性干擾。Oropeza等[33]提出多道奇異譜分析(Multichannel Singular Spectrum Analysis,MSSA)將地震數(shù)據(jù)變換到頻率域，然后對每個頻率切片構造塊Hankel矩陣，通過隨機奇異值分解對塊Hankel矩陣降秩，得到重建的地震數(shù)據(jù)。Kreimer等[34]將該框架直接運用于五維地震數(shù)據(jù)重建，通過高階奇異值分解(Higher Order Singular Value Decomposition)進行降秩。Abma等[35]提出凸集投影迭代方法，在頻率切片上進行傅里葉變換后，對幅值進行閾值處理，最后進行反傅里葉變換并插值。Jia等[36]提出正交矩陣匹配追蹤Hankel重建方法(Orthogonal Matrix Pursuit HankelReconstruction,OMPHR)，利用秩1正交匹配追蹤方法(Orthogonal Rank-one Matrix Pursuit)代替SVD分解進行降秩，大幅提升運算效率。Siahsar等[37]提出最優(yōu)奇異值收縮法，根據(jù)隨機矩陣理論產(chǎn)生最優(yōu)低秩估計量，并利用衰減因子約束奇異值，提高低秩估計的魯棒性[38]。Gao等[39]引入一種基于構建塊Toeplitz矩陣的快速降秩方法，運用Lanczos雙對角化算法代替截斷奇異值分解操作，提升了計算效率。Ma[40]提出基于構建紋理塊的低秩矩陣補全方法，通過三維紋理塊映射將三維數(shù)據(jù)變換到二維空間進行地震數(shù)據(jù)重建。

對于稀疏信號的恢復和低秩矩陣補全，一般建立秩最小化模型并將該模型凸松弛到核范數(shù)最小化模型。相對于凸核范數(shù)松弛模型，非凸Lp范數(shù)松弛模型可以以較低的采樣率恢復缺失數(shù)據(jù)，通過在迭代中設置權重約束奇異值，很好地保證了數(shù)據(jù)結構的低秩性[41-42]。Lu等[43]成功地將非凸Lp范數(shù)松弛模型應用到彩色圖像恢復上，并取得較好的效果。

基于降秩理論，本文提出非凸Lp范數(shù)Hankel重建(Non-convex LpNorm Hankel Reconstruction，NLPHR)方法對三維地震數(shù)據(jù)進行重建。非凸Lp范數(shù)比核范數(shù)更接近秩函數(shù)，因此通過求解最小Lp范數(shù)問題可以得到更接近于真實值的解。此外，MSSA方法和OMPHR方法均要求給定秩的大小，通過多次重復實驗選取合適的秩。本文方法不需要給定秩的大小，而是通過設置權重約束奇異值，迭代減小奇異值，保證了重建數(shù)據(jù)的低秩性。

1 基于非凸Lp范數(shù)的地震數(shù)據(jù)重建

1.1 三維地震數(shù)據(jù)低秩重建

考慮三維地震數(shù)據(jù)體D(t,x,y)∈RNt×Ny×Nx，其中Nt為時間采樣點數(shù)，Ny和Nx分別為x方向和y方向地震道數(shù)目。對每道數(shù)據(jù)進行傅里葉變換，得到頻率域數(shù)據(jù)D(ω,x,y)。取每一個頻率切片Dω∈RNy×Nx構成矩陣

(1)

首先對Dω的每行構造Hankel矩陣，以第j(1≤j≤Ny)行為例，格式如下

(2)

式中：Lx=floor(Nx/2)+1，其中floor(·)為向下取整函數(shù)；Kx=Nx-Lx+1。再對每一行的hj(1≤j≤Ny)進行Hankel矩陣化，建立塊Hankel矩陣

(3)

(4)

式中算子PH∶RNt×Ny×Nx→R(Lx×Ly)(Kx×Ky)表示將三維數(shù)據(jù)映射為二維塊Hankel矩陣。理論證明，同相軸數(shù)量無缺失的或無噪的地震數(shù)據(jù)是低秩的[33]，數(shù)據(jù)的缺失會增加矩陣的秩[40,44]。

1.2 NLPHR方法

低秩矩陣補全問題的一般形式為

(5)

(6)

對于地震數(shù)據(jù)重建問題，考慮建立非凸Lp范數(shù)松弛模型

(7)

假設1：h(x)∶R+→R+在[0,+∞)是凹函數(shù)且單調(diào)遞增，可以為非光滑函數(shù)；

假設2：d(x)∶Rm×n→R+是一個光滑函數(shù)，即其梯度為Lipschitz常數(shù)Lf，滿足

(8)

定義1：h(x)∶Rn→R是凹函數(shù)，如果對任意的b∈Rn，v滿足

h(b)≤h(a)+〈v,b-a〉

(9)

由非凸Lp范數(shù)h(x)=xp(0

(10)

根據(jù)假設1及定義1(式(9))，有

(11)

其中

(12)

由h(x)的增量性及超梯度的反單調(diào)特性，有

(13)

因為d(x)是一個Lf光滑函數(shù)，可將d(X)在X(k)處線性變換為

(14)

用式(11)右邊項代替式(7)的第一項，用式(14)右邊項代替式(7)第二項，得到式(7)解的更新值

(15)

求解式(15)等價于計算加權核范數(shù)的逼近算子，根據(jù)式(15)的非凸性及式(13),可得到式(15)的閉式解。根據(jù)

定理1[45]：對任意λ>0,Y∈Rm×n,0≤w1≤w2≤…≤wq,q=min(m,n)，Y=UΣVT是Y的奇異值分解，Sλw(Σ)=Diag[max(Σii-λwi,0)]。對于問題

(16)

其全局最優(yōu)解可通過加權奇異值閾值(WSVT)算法得到

X*=USλw(Σ)VT

(17)

通過求解式(15)更新X(k+1)后，根據(jù)

(18)

基于NLPHR的三維地震數(shù)據(jù)低秩重建流程(圖1)如下。

(1)首先將三維數(shù)據(jù)D(t,x,y)從時間域變換至頻率域，得到頻率域數(shù)據(jù)D(ω,x,y)，然后對每個頻率切片D(iΔω,x,y)分別構建相應的塊Hankel矩陣H。

(4)迭代停止后，從降秩得到的近似矩陣H*反對角線平均值獲取頻率域數(shù)據(jù)D*(ω,x,y)，將頻率域數(shù)據(jù)變換回時間域， 得到時間域數(shù)據(jù)D*(t,x,y)。

圖1 基于NLPHR的三維地震數(shù)據(jù)重建流程

2 實驗

通過對三維模擬地震數(shù)據(jù)和實際地震數(shù)據(jù)分別利用NLPHR法與MSSA方法、OMPHR方法進行對比。衡量重建質(zhì)量的指標為信噪比

(19)

2.1 模擬數(shù)據(jù)實驗

測試選用的模擬數(shù)據(jù)包含三條線性同相軸，x方向和y方向各有32道，道間距為1m。共有64個時間采樣點，采樣間隔為2ms。圖2a為原始模擬三維數(shù)據(jù)，其大小為64×32×32。每個塊Hankel矩陣大小為256×289。圖2b為隨機缺失40%道集的數(shù)據(jù)，圖2c為NLPHR重建結果。

為了更清楚地說明文中所提出方法的重建效果，分別取圖2中三維數(shù)據(jù)的二維剖面(y=19m)進行顯示(圖3)。圖3直觀地展示了NLPHR法重建前、后地震數(shù)據(jù)的吻合程度。圖4為OMPHR與MSSA方法重建效果對比，兩種方法中秩的大小均選取為3。NLPHR方法重建數(shù)據(jù)的信噪比為40.0dB，耗時3.8s； OMPHR和MSSA方法重建數(shù)據(jù)的信噪比分別為34.1dB和29.2dB，耗時分別為1.5s和94.5s。圖5為分別隨機缺失10%、20%、30%、40%、50%道集情況下NLPHR法、OMPHR法以及MSSA方法重建數(shù)據(jù)的信噪比曲線，可以看出，NLPHR方法重建信噪比明顯高于OMPHR法和MSSA法。

圖2 基于NLPHR法的三維模擬數(shù)據(jù)重建效果

圖3 基于NLPHR法顯示的重建剖面(y=19m)

2.2 實際數(shù)據(jù)實驗

實際三維疊后地震資料(圖6a)參數(shù)如下：x方向和y方向各有64道，道距為50m，時間采樣間隔為2ms，采樣點數(shù)為256。隨機缺失40%道集數(shù)據(jù)見圖6b，采用NLPHR法重建的地震數(shù)據(jù)如圖6c所示。圖7是從該三維數(shù)據(jù)體(圖6)中抽取的二維剖面(y=1000m)。從圖6和圖7的重構結果可知，基于NLPHR方法的三維數(shù)據(jù)重建方法能有效恢復缺失道集，較好地實現(xiàn)(含缺失道)地震數(shù)據(jù)的重建。圖8為隨機缺失40%道集情況下OMPHR法及MSSA法重構結果。OMPHR法和MSSA法的秩均取為20。由于數(shù)據(jù)體較大，實驗中實行分塊并行重建。將大小為256×64×64的數(shù)據(jù)均分成四塊，每塊大小為256×32×32，對四塊同時進行重建。NLPHR法重建數(shù)據(jù)的信噪比為11.9dB，耗時81.7s；OMPHR法重建數(shù)據(jù)的信噪比為9.3dB，耗時66.6s；MSSA重建數(shù)據(jù)的信噪比為7.7dB，耗時203.2s?？梢钥闯?，NLPHR法較好地恢復了缺失道集，而OMPHR法和MSSA法均未能完全恢復缺失道集。圖9為不同缺失率下三種方法重建數(shù)據(jù)的信噪比對比?？梢钥闯?，NLPHR法具有較高的輸出信噪比，比OMPHR法和MSSA法更具有穩(wěn)健性。表1為三種數(shù)據(jù)重建方法在不同缺失率下三維疊后地震資料上重建數(shù)據(jù)的信噪比及計算耗時對比。從表1可以看出，NLPHR法重建數(shù)據(jù)的信噪比高于OMPHR法和MSSA法，耗時比MSSA方法短。圖10和圖11分別為缺失40%數(shù)據(jù)下該三維疊后數(shù)據(jù)體在參數(shù)p和η不同取值時重建數(shù)據(jù)信噪比及耗時對比。由圖可見，參數(shù)p在(0.5,1)范圍內(nèi)重建數(shù)據(jù)的信噪比較高；p=0.6時重建信噪比最高，且重建耗時相對較短。由圖11可以看出，η在(0.1,0.8)范圍內(nèi)，重建信噪比隨η呈現(xiàn)增大趨勢；在(0.8,0.9)范圍內(nèi)呈下降趨勢。圖12為數(shù)據(jù)缺失40%時三維疊后數(shù)據(jù)在三種重建方法下的SNR收斂曲線及收斂時間。由圖可知，三種方法均較快收斂，其中NLPHR法收斂最快，重建數(shù)據(jù)的信噪比高于OMPHR法和MSSA法，迭代時間比MSSA法短。

圖4 重建剖面(y=19m)結果對比

圖5 三維模擬數(shù)據(jù)在不同缺失率下不同方法重建數(shù)據(jù)信噪比對比

圖6 基于NLPHR法的三維疊后地震數(shù)據(jù)重建結果

圖7 基于NLPHR法的重建結果二維剖面(y=1000m)顯示

圖8 疊后數(shù)據(jù)不同方法重建結果的二維剖面(y=1000m)對比

圖9 三維疊后數(shù)據(jù)在不同缺失率下、不同方法重建數(shù)據(jù)信噪比對比

表1 三維疊后數(shù)據(jù)在不同缺失率下三種方法重建數(shù)據(jù)信噪比及耗時對比

圖10 隨機缺失40%道集的疊后數(shù)據(jù)在參數(shù)p取不同值下重建數(shù)據(jù)信噪比(左)及計算耗時時間(右)曲線

圖11 隨機缺失40%道集的疊后數(shù)據(jù)在參數(shù)η取不同值下重建數(shù)據(jù)信噪比(左)及重建時間(右)曲線

圖12 隨機缺失40%道集的疊后數(shù)據(jù)重建數(shù)據(jù)的SNR (左)及收斂時間(右)曲線

3 結束語

MSSA法通過隨機奇異值分解直接降秩，得到的解往往不是最優(yōu)解，并且在低采樣率下重建效果欠佳。針對此問題，文章提出非凸Lp范數(shù)Hankel重建方法(NLPHR)，對三維地震數(shù)據(jù)進行重建。該方法將地震數(shù)據(jù)每個頻率切片轉(zhuǎn)換成塊Hankel矩陣，再運用NLPHR法對塊Hankel矩陣進行降秩處理，從而實現(xiàn)地震數(shù)據(jù)低秩重建。三維模擬地震數(shù)據(jù)和實際地震數(shù)據(jù)實驗證明，與MSSA法和OMPHR法相比，本文方法數(shù)據(jù)重建信噪比最高，耗時較少，效果最優(yōu)。