非凸Lp范數三維地震數據重建

劉 群 付麗華 張婉娟

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

0 引言

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

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

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

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

1 基于非凸Lp范數的地震數據重建

1.1 三維地震數據低秩重建

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

(1)

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

(2)

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

(3)

(4)

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

1.2 NLPHR方法

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

(5)

(6)

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

(7)

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

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

(8)

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

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

(9)

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

(10)

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

(11)

其中

(12)

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

(13)

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

(14)

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

(15)

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

定理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)

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

X*=USλw(Σ)VT

(17)

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

(18)

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

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

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

圖1 基于NLPHR的三維地震數據重建流程

2 實驗

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

(19)

2.1 模擬數據實驗

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

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

圖2 基于NLPHR法的三維模擬數據重建效果

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

2.2 實際數據實驗

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

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

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

圖6 基于NLPHR法的三維疊后地震數據重建結果

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

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

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

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

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

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

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

3 結束語

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