基于壓縮感知的L1范數(shù)譜投影梯度算法地震數(shù)據(jù)重建

蘭天維,韓立國,張 良

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

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

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

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

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

1 方法理論

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

1.1 壓縮感知理論

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

(1)

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

(2)

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

(3)

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

(4)

該問題的求解方法為:

(5)

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

(6)

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

1.2 contourlet變換

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

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

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

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

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

圖1 contourlet變換流程

1.3 SPGL1算法

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

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

(7)

式中:正參數(shù)σ為對(duì)噪聲水平的估計(jì),σ=0時(shí)為BP問題的解。

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

(8)

式中:τ為閾值。

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

(9)

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

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

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

(10)

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

(11)

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

(12)

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

SPGL1算法流程如下:

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

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

5) 返回步驟2)。

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

2 模型數(shù)據(jù)算例

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

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

(13)

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

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

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

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

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

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

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

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

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

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

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

表2 3種重建算法重構(gòu)合成地震記錄的信噪比

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

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

圖7 計(jì)算時(shí)間隨迭代次數(shù)變化的曲線

3 實(shí)際地震數(shù)據(jù)處理及分析

選擇復(fù)雜地區(qū)實(shí)際地震數(shù)據(jù)進(jìn)行處理和分析。

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

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

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

圖9 在contourlet稀疏變換基下3種算法對(duì)實(shí)際地震數(shù)據(jù)重建的結(jié)果對(duì)比a OMP算法; b GPSR算法; c SPGL1算法

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

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

表3 三種重建算法重構(gòu)實(shí)際地震記錄的信噪比

4 結(jié)束語

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

合成地震數(shù)據(jù)和實(shí)際地震數(shù)據(jù)測試結(jié)果驗(yàn)證了該方法的可行性和有效性。在相同重建算法下,相比傅里葉變換與小波變換,在高精度重建能量較強(qiáng)的同相軸之余,contourlet變換清晰刻畫了由復(fù)雜地質(zhì)構(gòu)造產(chǎn)生的能量較弱的反射波、繞射波信息;在相同稀疏變換基下,SPGL1算法實(shí)現(xiàn)了較高精度的重建,且計(jì)算效率最高。在實(shí)際地震數(shù)據(jù)測試中,相比OMP算法與GPSR算法,SPGL1算法重建精度高、速度快、殘差小,同時(shí)兼顧了對(duì)能量較弱的反射波、繞射波的重建。

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