蘭天維,韓立國,張 良

(吉林大學地球探測科學與技術學院,吉林長春130026)

隨著我國油氣勘探程度的加深,勘探區(qū)域的范圍與地質構造的復雜程度愈發(fā)加大,采集到的地震數據的規(guī)模與復雜度不斷增加。復雜、大范圍的勘探區(qū)域與地形、障礙物等環(huán)境因素以及儀器等經濟因素一起,加劇了地震數據的不規(guī)則性或稀疏分布,這將不同程度地影響地震數據的后續(xù)處理效果,因此需要對不完整采樣下的缺失地震道進行插值處理,即地震數據重建。目前地震數據重建的方法主要分為3類:基于濾波器的策略、基于波動方程算子的方法以及基于某種變換的重建技術?；跒V波器的重建技術利用褶積插值濾波器將不規(guī)則數據當做規(guī)則數據處理,常用的方法是預測誤差濾波法[1-2],此類方法會造成較大誤差,其結果也具有不確定性?；诓▌臃匠痰牟逯捣椒ㄖ饕ㄟ^正演與反演算子來迭代求解一個反問題,如NMO、反DMO、方位角時差校正AMO等[3-4],此類方法需要地下結構的先驗信息,計算量大,對采樣率要求較高。基于某種變換的方法通常是對地震數據進行某種變換,然后在變換域對地震數據進行重建。主要有傅里葉變換、小波變換[5]、Randon變換[6]、curvelet變換[7]、Seislet變換、contourlet變換[8]等方法。

壓縮感知理論[9-10]是一種以低于Nyquist采樣率采樣的理論,其本質是當信號具有稀疏性或者在某個變換域可以稀疏表示時,利用一個與稀疏變換基不相關的觀測矩陣,將該信號從高維映射到低維,然后利用稀疏促進算法求解最優(yōu)化問題,實現數據的高概率重構。地震信號可以通過某種變換進行稀疏表示滿足了壓縮感知理論應用的條件?；趬嚎s感知的地震數據重建常用的變換方法主要有離散余弦變換(DCT)、傅里葉變換、小波變換、curvelet變換和學習型超完備字典等。DCT與傅里葉變換是全局變換,無法有效識別地震信號局部時頻特征;短時傅里葉變換時頻局部化能力有限,無法滿足不同頻率與時間的信號分辨要求;小波變換雖然彌補了短時傅里葉變換的缺點,但是無法很好地描述二維信號中的線性特征;curvelet變換具有方向識別能力,可近似最優(yōu)地表示二維信號中的曲線特征,然而其冗余度較高,計算效率低;學習型超完備冗余字典的變換能夠很好地稀疏表示信號,但迭代次數多,訓練時間長。相較于傅里葉變換與小波變換,contourlet變換以其方向性、局部性與各向異性,可較好地處理信號的線性、曲線特征,善于處理信號的奇異性,捕捉地震記錄中由直達波、反射波、繞射波等波場構成的輪廓信息,從而能夠較清晰地描述信號中因復雜構造產生的各種波場。與學習型超完備字典相比,contourlet變換運算時間短,算法簡單,更加符合實際應用中的計算效率要求。

在基于壓縮感知的地震數據重建方法中,重建算法是重要的一環(huán)。近年來基于壓縮感知的重建算法主要有3種,分別是貪婪算法、組合算法和凸優(yōu)化算法。其中貪婪算法通過迭代尋找一組與觀測值匹配的最稀疏原子,實現信號重建,包括匹配追蹤(MP)[11]和正交匹配追蹤(OMP)[12-13]等。組合算法要求原始信號能夠支持快速分組測試重建,包括組測試和數據流草圖[14-15],該算法需要對觀測矩陣進行進一步設計。凸優(yōu)化算法是將非凸問題轉化為凸問題求解,尋找信號的逼近,該算法的一個極為突出的優(yōu)點是,在目標函數為嚴格凸函數時,用此算法求解得到的局部極大值(極小值)就是全局最大值(最小值),并且所求結果只有一個。凸優(yōu)化算法包括基追蹤[16]、迭代閾值[17]、內點法和梯度投影法[18-19]。其中基追蹤算法為凸優(yōu)化算法中的基本算法,由于計算量大,基本上用于一維信號處理;迭代閾值法運算簡單,但收斂速度慢;內點法一般用于中小規(guī)模問題;投影梯度算法[20-22]應用簡單并且適合大規(guī)模與復數域問題,其中SPGL1收斂速度較普通的梯度投影法快[19]。

綜合以上3種重建方法,考慮各稀疏變換與重建算法的精度及計算效率,本文提出了一種基于壓縮感知的地震數據重建方法,該方法以contourlet變換為稀疏變換基,以L1范數譜投影梯度法為重建算法。該方法應用于模型數據和實際地震數據的處理結果驗證了其在大規(guī)模、復雜地震數據情況下的有效性。

1 方法理論

基于壓縮感知的L1范數譜投影梯度算法地震數據重建方法由三部分組成:壓縮感知理論、contourlet變換與SPGL1算法,其中壓縮感知理論為主體,contourlet變換為壓縮感知的稀疏變換階段,SPGL1為壓縮感知的重建階段。

1.1 壓縮感知理論

壓縮感知理論為地震數據的壓縮采樣提供了可能。在該理論框架下,地震數據的采樣表達式為:

(1)

式中:f∈Rn為原始地震數據;y∈Rm(m?n)為采樣后地震數據;φ∈Rm×n為采樣矩陣。地震數據的重建即將y恢復為f。若f在某個變換域中可稀疏表示,即:

(2)

式中:φ為稀疏變換算子;x為f的稀疏表示系數。將(1)式和(2)式結合,那么壓縮感知理論框架下的數據重建可以表述成:

(3)

式中:φH為稀疏逆變換算子。以感知矩陣A代替φφH,則(3)式可寫為:

(4)

該問題的求解方法為:

(5)

此方法也稱之為L0優(yōu)化問題,但是由于離散性與L0范數的非凸性,它是一個NP難題。為了解決此問題,人們用L1替代L0[23]:

(6)

在壓縮感知中,為了保證(4)式的可解性,對矩陣A進行了限定,即矩陣A應滿足有限等距性(restricted lsometry property,RIP)[23]。矩陣A是否滿足RIP,由φ與φH是否相干決定,兩者不相干時,A大概率滿足RIP。在地震勘探中,常常把缺失的地震道看做0,未缺失的地道看做1[24],由此組成的對角矩陣被作為采樣矩陣φ時,φ與φH不相干,滿足了矩陣A的限定條件。因此,本文將此對角矩陣作為采樣矩陣。

1.2 contourlet變換

隨著地震勘探的深入,采集到的信號中包括越來越多的反射波、繞射波等輪廓信息,但是傅里葉變換與小波變換無法稀疏表示輪廓線,難以適應包含較多輪廓信息的復雜地震信號,因此,DO[25]提出了contourlet變換,這是一種二維可采集信號內在幾何結構的變換,具有靈活的方向性和各向異性。它用類似于輪廓段的基結構來逼近幾何結構,基的支撐區(qū)間是長寬比隨尺度變化的“長條形”結構,可以對分段函數實現最優(yōu)近似。

在壓縮感知理論下,稀疏表示地震信號的contourlet變換步驟如下:

1) 使用拉普拉斯金字塔濾波器(LP)對地震信號進行多尺度分解,捕獲不同尺度上的邊緣奇異點。每一級LP分解首先對上一級低頻信息進行下采樣產生本級低頻信息,然后對此低頻信息進行上采樣并與上一級低頻信息相減得到高頻信息,最后對低頻信息進行LP分解迭代即實現了多尺度分解。

2) 對由LP分解產生的高頻信息使用方向濾波器組(DFB)進行方向分解,并連接同方向的奇異點,得到contourlet變換系數。如DO[25]提出的contourlet變換中的DFB,將雙通道梅花濾波器組與剪切算子結合,采用簡化的二叉樹分解方法,可獲得較好的方向信息。

圖1為contourlet變換流程圖[26]。

圖1 contourlet變換流程

1.3 SPGL1算法

針對凸優(yōu)化問題一般具有規(guī)模大、目標函數非光滑的特點,本文采用L1范數譜投影梯度算法(SPGL1)來解決凸問題,得到重建的信號。

(6)式又被稱為基追蹤問題(basis pursuit,BP),其中A為m×n的矩陣,y為m列的列向量。當數據含有噪聲或數據不完整時,基追蹤問題變?yōu)榛粉櫲ピ?basis pursuit denoise,BPDN):

(7)

式中:正參數σ為對噪聲水平的估計,σ=0時為BP問題的解。

BPDN僅為L1范數正則化最小二乘問題中的一種,它還包括Lasso問題,其表達式如下:

(8)

式中:τ為閾值。

BPDN常用在懲罰最小二乘問題上,其形式是:

(9)

其中,λ是一個與Lasso問題中的約束拉格朗日乘數和BPDN問題中的約束乘數的倒數有關的參數。參數σ,λ,τ的適當選擇,使得BPσ問題、LSτ問題、QPλ問題的解在某些情況下一致。在σ已知的情況下,λ可用同倫算法求取[19]。

因此,SPGL1算法的思想是將BPDN問題轉化為一系列的Lasso子問題,通過譜投影梯度法(SPG)求解Lasso問題,以得到BPDN問題的解。

令xτ為LSτ問題的最優(yōu)解,對于任意τ≥0,Lasso問題的最優(yōu)值為:

(10)

Pareto曲線定義了殘差r的L2范數與解x的L1范數之間的平衡,BPDN與Lasso是同一條曲線的兩個不同特征,可以用參數τ參數化Pareto曲線。對于定義域內所有的點,ψ都是連續(xù)可微的,這樣可以使用牛頓法求解(11)式:

(11)

(11)式定義了一系列正則化參數τk→τσ,可表示為:

(12)

式中:ψ′(τk)為ψ(τk)的導數。這樣Lasso問題的相關解xτσ收斂到xσ,參數τσ使得BPDN與Lasso得到相同解。

SPGL1算法流程如下:

1) 初始化參數,步長α0∈[αmin,αmax],充分下降常數γ∈(0,1),初始迭代x0←Pτ[x],初始殘差r0←y-Ax0,初始梯度方向g0←-ATr0,l←0,整數線性搜索步長M≥1。

2) 計算對偶間隙δl←‖rl‖2-(yTrl-τ‖gl‖∞)/‖rl‖2,收斂則退出,轉步驟6)。

5) 返回步驟2)。

6) 輸出xτ←xl,rτ←rl。

2 模型數據算例

2.1 相同重建算法下合成地震記錄測試及分析

本文采用信噪比來衡量地震數據的重建效果:

(13)

選擇大規(guī)模復雜合成地震記錄測試本文方法的重建效果,速度模型(Marmousi模型)如圖2所示。

圖2 大規(guī)模復雜合成數據的速度模型(Marmousi模型)

該模型包括1043×302個網格點,空間網格間距為10m,模型上部與2500m深度處共有4個斷層,斷層的存在使得地震記錄有多處繞射波,增大了地震數據重建的難度。

在地表激發(fā)地表接收觀測系統(tǒng)下,以主頻為20Hz的Ricker子波為震源,在該模型地表4000m處激發(fā),放置1043個檢波器,時間采樣間隔為1ms,得到原始地震記錄如圖3a所示。

在壓縮感知理論下采用同一種重建算法,不同的稀疏變換基對隨機缺失50%的地震記錄(圖3b)進行重建。圖4a至圖4c分別為以傅里葉變換、小波變換、contourlet變換為稀疏變換基,應用SPGL1算法進行數據重建的結果,參數σ=0.05,最大迭代次數為30[27]。表1對比了3種稀疏變換基采用相同重建算法的信噪比(RSN)。由圖3可以看到,合成地震數據中含有較多直達波、反射波、繞射波等波場輪廓信息。由圖4可見,由于傅里葉變換識別局部時頻特征的能力有限,無法有效處理地震記錄中的輪廓信息,在稀疏變換時只能勉強識別能量較強的同相軸,無法較好識別震源處地震記錄特征,重建結果存在噪聲,能量較弱的反射波、繞射波等難以識別(圖4a中紅色箭頭處);小波變換相較于傅里葉變換可以較為清晰地識別同相軸,然而依舊無法較好地刻畫各種輪廓信息,因此重建結果仍然有噪聲存在,能量較弱的反射波、繞射波無法清晰刻畫(圖4b);而采用contourlet變換得到的重建地震記錄質量好,同相軸連貫,反射波、繞射波均清晰可見(圖4c)。

圖3 原始地震記錄(a)與隨機缺失50%的地震記錄(b)

圖4 不同稀疏變換基下應用SPGL1算法的重建效果a 傅里葉變換; b 小波變換; c contourlet變換

表1 3種稀疏變換基重構合成地震記錄的信噪比

2.2 相同稀疏變換基下合成地震記錄測試及分析

采用基于壓縮感知和L1范數譜投影梯度算法重建地震數據的方法對圖3b地震數據進行重建。其中SPGL1算法參數選擇為σ=0.01,迭代次數為50次[27]。圖5對比了contourlet稀疏變換基下OMP、GPSR、SPGL1算法重建的合成地震記錄。圖6 為相應的殘差,表2給出了3種重建算法重構合成地震記錄的信噪比。由圖5、圖6和表2可以看到,3種算法均可較好地完成合成地震記錄的重建,而且差別不大。但是在相同稀疏變換基下,采用OMP算法進行重建的50次平均CPU時間為8164s,GPSR算法為14532s,而SPGL1算法為991s,遠遠小于前兩種算法重建所用的時間。圖7為CPU時間隨迭代次數變化的曲線,可見,隨著迭代次數的增加,SPGL1算法不僅用時較少,而且整體變化不大;GPSR算法用時漸長;OMP算法隨著迭代次數的增加,CPU時間大幅增加。

注意到在相同稀疏變換基下對合成數據進行重建時,雖然3種算法重建結果差別不大,但是OMP算法信噪比低于其它兩種算法,這是由于OMP算法通過迭代得到的局部最優(yōu)解并不一定是全局最優(yōu)解,凸優(yōu)化算法的特點是利用L1范數求解優(yōu)化問題時得到的局部最優(yōu)解為全局最優(yōu)解。GPSR與SPGL1算法重建信噪比相同。

表2 3種重建算法重構合成地震記錄的信噪比

圖5 3種算法重建的合成地震記錄對比a OMP算法; b GPSR算法; c SPGL1算法

圖6 3種重建算法重建結果的殘差a OMP算法; b GPSR算法; c SPGL1算法

圖7 計算時間隨迭代次數變化的曲線

3 實際地震數據處理及分析

選擇復雜地區(qū)實際地震數據進行處理和分析。

圖8a為原始記錄,共228道,每道3001個檢波器,道間距為25m,時間采樣間隔為2ms。SPGL1算法參數與2.2節(jié)中SPGL1算法參數的選擇相同。

針對原始記錄隨機缺失50%的記錄(圖8b),使用基于壓縮感知和L1范數譜投影梯度算法重建地震數據的方法進行地震數據重建。圖9為在contourlet稀疏變換基下,采用OMP、GPSR、SPGL1算法的重建結果,圖10為相應的殘差,表3 為3種算法重建實際地震數據的信噪比。可以看出,波場中直達波、繞射波發(fā)育,3種算法均可實現地震記錄的重建,其中OMP算法重建效果與SPGL1算法類似。但對比相應的殘差(圖10a和圖10c)可知,OMP算法重建信噪比較低,留有較多殘差,其運行50次的平均CUP時間為39365s,而同樣情況下SPGL1算法運行50次的平均CPU時間為476s。GPSR算法的重建結果信噪比低,只能較好地識別能量較強的同相軸,能量較弱的反射波基本為噪聲所覆蓋(圖9b紅色箭頭處),其重建用時達10506s,遠不能滿足實際工作中對于精度與效率的要求。

圖8 實際地震記錄(a)與隨機缺失50%的地震記錄(b)

圖9 在contourlet稀疏變換基下3種算法對實際地震數據重建的結果對比a OMP算法; b GPSR算法; c SPGL1算法

圖10 3種重建算法的殘差對比a OMP算法; b GPSR算法; c SPGL1算法

對表3進行分析并與表2比較發(fā)現,OMP的重建精度低于SPGL1算法,但卻高于GPSR算法。這是由于實際數據中含有大量噪聲,在對此數據進行重建時,GPSR算法相對于SPGL1算法的魯棒性較差,受噪聲影響大。在對實際數據進行濾波后,再采用3種算法進行重建后發(fā)現,OMP算法重建精度依然較低,信噪比為23.8,GPSR與SPGL1算法的信噪比相差不大,分別為26.2與26.8,與合成數據重建結果較為一致。

表3 三種重建算法重構實際地震記錄的信噪比

4 結束語

本文基于壓縮感知理論,以contourlet變換做為稀疏變換基,結合SPGL1重建算法,提出了一種基于壓縮感知的L1范數譜投影梯度算法地震數據重建方法。

合成地震數據和實際地震數據測試結果驗證了該方法的可行性和有效性。在相同重建算法下,相比傅里葉變換與小波變換,在高精度重建能量較強的同相軸之余,contourlet變換清晰刻畫了由復雜地質構造產生的能量較弱的反射波、繞射波信息;在相同稀疏變換基下,SPGL1算法實現了較高精度的重建,且計算效率最高。在實際地震數據測試中,相比OMP算法與GPSR算法,SPGL1算法重建精度高、速度快、殘差小,同時兼顧了對能量較弱的反射波、繞射波的重建。

對于大規(guī)模復雜數據的重建,精度與效率是兩個需要考量的重要因素。近年來學者們在壓縮感知理論框架下提出了許多卓有成效的地震數據重建方法,然而面對地質情況復雜、勘探面積大、地震數據量劇增的情況,今后仍需要在進一步提高重建精度的同時兼顧計算效率方面展開相關研究。