石 敏 朱震東 路 昊 朱登明 周 軍

(*華北電力大學控制與計算機工程學院 北京102206)

(**中國科學院計算所前瞻研究實驗室 北京100190)

(***中國石油集團測井有限公司 西安710077)

0 引言

地震勘探是油氣勘探的重要手段。隨著油氣勘探的不斷深入,勘探目標區(qū)域的結構和周圍環(huán)境越來越復雜,使得收集的數(shù)據(jù)不規(guī)則甚至不完整,從而對后續(xù)分析解釋及油氣判斷產生了影響。同時,這也意味著更高的勘探成本。如果能夠基于低密度采樣數(shù)據(jù)重建出準確的高密度地震數(shù)據(jù),便能夠改善現(xiàn)有勘探方法,有效地利用獲得的數(shù)據(jù)集可為估算地下石油儲量的形成提供更可靠的支持,并從一定程度上降低勘探成本。地震數(shù)據(jù)在變換域中呈現(xiàn)稀疏分布,并且在時域和空域上也顯示出很強的相關性,這為地震數(shù)據(jù)重建提供了可能。

在這項工作中,針對不完整地震數(shù)據(jù)的壓縮感知重建方法,本文分析了地震勘探和測井勘探的數(shù)據(jù)特點,研究了一種時空約束模型,對傳統(tǒng)壓縮感知模型進行了改進。本文通過內核奇異值分解(kernel singular value decomposition,K-SVD)算法來訓練超完備字典,并使用改進的稀疏自適應匹配追蹤算法(sparsity of adaptive matching pursuit,SAMP)解決相應的優(yōu)化問題,從而完成壓縮感知重建。最終,進行了大量的對比實驗,驗證了本文算法的重建效果及效率。

總而言之,本文的貢獻如下。

(1) 分析了地震數(shù)據(jù)的特征,改進了經(jīng)典的壓縮感知模型,添加時空相關信息作為壓縮感知模型的約束。

(2) 改進了稀疏自適應匹配追蹤算法,增加了初始稀疏性估計和可變步長的策略,確保了重建精度的同時,提高了算法效率。

(3) 在真實地震數(shù)據(jù)和微電阻率成像數(shù)據(jù)上實現(xiàn)了本文算法,證明了其出色的重建能力和泛化能力。

1 相關工作

基于地震數(shù)據(jù)的重建方法,通常被劃分為以下3 種:第一種方法是基于預測濾波的,通常通過高斯窗口對不規(guī)則樣本數(shù)據(jù)進行插值,這種方法導致很多錯誤出現(xiàn)[1-2]。第二種方法基于波動方程,通過正反演算子解決一個反問題。這種方法結合地下結構的先驗信息,并依據(jù)波傳播的物理特性實現(xiàn)重建,這種方法通常計算量很大[3]。第三種是基于變換的方法,首先對地震數(shù)據(jù)進行變換,進而在變換域中通過迭代求解等方法實現(xiàn)地震數(shù)據(jù)的重建。由于具備穩(wěn)定性和可解釋性,該方法得到了較為廣泛的發(fā)展和應用[4-5]。Abma 等人[6]在地震重建領域引入了壓縮感知算法[7-8]。Wen 等人[9]指出,在基于壓縮感知方法的重建中,影響地震數(shù)據(jù)重建效果及效率的3 個主要因素是稀疏變換,迭代算法以及閾值模型。Bora 等人[10]提出了一種基于生成網(wǎng)絡的壓縮感知算法,該算法實現(xiàn)了基于更少的數(shù)據(jù)進行重建,但是帶來了諸如訓練困難和精度不足的問題。

稀疏變換對重建效果有很大的影響,其包括離散余弦變換(DCT),傅立葉變換(Fourier)和超完備字典[11]等。在壓縮感知的應用領域中,DCT 和Fourier 較為常見,但這兩種變換都難以有效地識別局部特征。短時傅立葉變換改善了這一問題,但是對于諸如地震數(shù)據(jù)之類的復雜信號,在不同時間的波形變化很大,短時傅立葉變換的時頻局部化能力仍然有限。小波分析采用并完善了這種局部化思想[12],但是不能識別方向。直到后來,提出了曲波(Curvelet)變換[13-14],在保證多尺度識別能力情況下優(yōu)化了多方向識別能力。因此在地震數(shù)據(jù)重建領域,曲波變換作為最佳稀疏變換方法之一而被經(jīng)常應用[15]。剪切波變換[16]對多方向識別能力進一步發(fā)展,可以使地震信號的稀疏表示更加稀疏。

目前,基于上述變化的壓縮感知理論已被用于地震數(shù)據(jù)重建[17-19],取得了良好的數(shù)據(jù)重建效果。但是,上述變換方法有一個共同點,變換基都是預先選擇的,并不一定能夠適合特定場景數(shù)據(jù)本身的特征。由于實際的地震數(shù)據(jù)一般非常復雜,通常會涉及多種類型的元素,預設的稀疏基很難對復雜的地震數(shù)據(jù)進行有效變換。在相關研究中,使用內核奇異值分解算法來學習超完備字典,能夠結合場景數(shù)據(jù),訓練得到更適應數(shù)據(jù)特征的字典作為稀疏變換基[20-22]。在最近的研究中,基于字典學習的地震數(shù)據(jù)處理展示了比傳統(tǒng)字典更強的稀疏重建能力[23-25]。

迭代算法不僅會影響到整體的重建效果,還極大程度上決定了算法的效率。由于早期迭代算法(如內點算法[26]和梯度投影法[27])的缺點,提出了重建貪婪算法,包括傳統(tǒng)的貪婪算法,例如匹配追蹤算法和正交匹配追蹤算法(OMP)[28-29],以及基于傳統(tǒng)算法改進得到的分段正交匹配追蹤[30]和規(guī)則化正交匹配追蹤[31]等,但是這些方法通常只應用于信號具有較低的稀疏度的場景。子空間追蹤算法(SP)[32]和壓縮采樣匹配追蹤算法(CoSaMP)[33]是兩種具有良好性能的相似算法。Blumensath[34]提出了迭代硬閾值(IHT)及其改進算法,進一步提升了重建效果。但是,這些傳統(tǒng)的貪婪算法都需要預先獲得信號稀疏性。對此,Thong 等人[35]提出了一種自適應估計信號稀疏性的迭代算法,即稀疏自適應匹配追蹤(SAMP)算法。解決了估計信號稀疏性的問題,但是并沒有得到最佳的重建精度,并且算法效率較低。

綜上所述,研究基于壓縮感知的地震數(shù)據(jù)重建,主要在于研究如何基于地震數(shù)據(jù)的相關特征,構造稀疏效果更好的變換基以及該場景下更高效快速的迭代算法。另外,傳統(tǒng)的壓縮感知框架下的重建是針對單幀數(shù)據(jù)進行的,而地震數(shù)據(jù)在時域和空域上具有相關性,在重建地震數(shù)據(jù)時,僅考慮單幀數(shù)據(jù)的信息,而忽略了連續(xù)幀之間的相關信息,這將會從一定程度上影響地震數(shù)據(jù)的重建效果。因此本文也對此也進行了研究。

2 基于時空約束的壓縮感知

2.1 地震數(shù)據(jù)的時空相關性

地震信號是復雜的,甚至瞬態(tài)特征也是不穩(wěn)定的[36]。由于巖層的密度差異的存在,在不同巖層的交界處會存在地震波反射和地震波折射的現(xiàn)象,因此地震波數(shù)據(jù)與時間域相關,不同深度的地質信息可以由探測器在不同時間點獲取地震波來表示。而震波在同一層的巖石中連續(xù)傳遞,因此地震波數(shù)據(jù)也與空間域有關,不同空間位置的地質信息可以由波器延測線的水平分布來表示。

目前,傳統(tǒng)的壓縮感知算法能夠基于不完整的地震數(shù)據(jù)ym∈RM,通過重建得到完整的地震數(shù)據(jù)fm∈RN(M ＜N)。然而,傳統(tǒng)的壓縮感知算法僅滿足了單幀數(shù)據(jù)重建的合理性,而忽略了地震數(shù)據(jù)相鄰幀之間的相關性。如圖1 所示,利用傳統(tǒng)的壓縮感知算法,基于單幀信息進行數(shù)據(jù)重建,重建結果中的高頻區(qū)域具有較為明顯的“失真”現(xiàn)象。因此,在利用壓縮感知算法對地震數(shù)據(jù)進行重建時,有必要利用數(shù)據(jù)相鄰幀之間的相關信息對重建進行優(yōu)化或約束。

圖1 單幀數(shù)據(jù)重建結果

2.2 基于時空約束的壓縮感知算法

壓縮感知算法使用觀測矩陣用于描述采樣數(shù)據(jù)和完整地震數(shù)據(jù)之間的關系。

其中,Φm∈RM×N是觀測矩陣,而m是地震道號。利用稀疏變換基φ進行稀疏,完整的地震數(shù)據(jù)fm可以進而表示為

其中,xm是稀疏向量,將兩式結合在一起,可以得到壓縮感知的基本表示形式:

在傳統(tǒng)的壓縮感知算法基礎上進行了修改,增加時空相關性約束,使其能夠更好地重建連續(xù)多幀地震數(shù)據(jù)。設{fm-n,fm-n+1,fm-n+2,…,fm} 為重建的連續(xù)n+1 個數(shù)據(jù)幀,并定義R為連續(xù)幀之間變化的能量損失。

可以看出,R可以衡量幀之間的差異,并與相關性呈負相關。將式(2)帶入式(6),可以得到:

目標是找到R的最小值,這等同于求解目標。對于稀疏數(shù)據(jù),最小L1 范式和最小L0 范式從一定程度上是等效的。因此,可以將目標轉換為最小L0 范式的解。

從而可以為式(4)添加時空相關信息R,使重建的數(shù)據(jù)具有盡可能小的變化損失:

其中,λ代表時空相關損失的比例。還可以在式(9)上進行等效轉換,以使加入時空信息的模型與壓縮感知的表示相同,便于使用迭代算法來重建數(shù)據(jù)。本文將λ設為1,然后進行等效轉換得到:

3 地震數(shù)據(jù)重建

本節(jié)對地震數(shù)據(jù)重建算法進行實現(xiàn)。地震數(shù)據(jù)的特征比較復雜,不同地質環(huán)境下的地震數(shù)據(jù)也有很大差異。非自適應變換基不能很好地表示地震特征。因此,選擇利用K-SVD 算法來訓練超完備字典,以此用重建所需的稀疏變換基。本文還通過添加初始稀疏性估計和變步長更新策略來改進SAMP算法?；诘? 節(jié)中的優(yōu)化模型,使用改進的SAMP 算法來求解稀疏矩陣。將稀疏矩陣和超完備字典相乘以獲得重建數(shù)據(jù)。完整的流程圖如圖2 所示。

圖2 地震數(shù)據(jù)重建流程圖

3.1 基于K-SVD 的超完備詞典學習

K-SVD 算法的本質是迭代思想,每次迭代使用范數(shù)稀疏約束跟蹤計算稀疏系數(shù),并利用奇異值分解算法更新字典原子,使得稀疏系數(shù)與字典能夠得到同步更新,最終可以根據(jù)稀疏約束來自適應地訓練出超完備字典。

假設,矩陣Y=是N個完整的地震數(shù)據(jù)切片的集合,代表了字典學習過程中的訓練數(shù)據(jù),矩陣X=是N個與訓練數(shù)據(jù)Y=所對應的稀疏表示系數(shù)向量集合,矩陣φ∈Rn×K表示超完備字典。那么超完備字典訓練的過程就可以表示為一個優(yōu)化問題:

其中,T0是稀疏表示系數(shù)中非零元素的最大數(shù)量。超完備字典訓練的整體步驟如下。

步驟1利用某種變換基進行字典初始化。

步驟2根據(jù)已知字典φ,使用迭代算法求解每個樣本yi的稀疏系數(shù)向量xi,即

步驟3更新字典φ。如果滿足收斂條件抑或達到預先設置的迭代次數(shù)上限,則獲得最終字典φ,否則轉向步驟2。設向量dk為要更新的字典φ的第k列原子,此時樣本集的分解形式可表示為

其中,向量是與dk相對應的X中第k行向量;矩陣Ek代表提取dk之后的誤差矩陣。

3.2 基于自適應動態(tài)步長SAMP 算法的重建

在地震數(shù)據(jù)重建領域通常情況下待重建數(shù)據(jù)的稀疏性是未知的,而在這種情況下,SAMP 算法可以通過設置固定步長s動態(tài)更新稀疏度來達到最優(yōu)。通常設定稀疏度為最小值,并在每次迭代中,根據(jù)殘值判斷是否需要增加稀疏度,逐漸逼近信號實際的稀疏度,從而得到稀疏度的最佳估計,最終實現(xiàn)重建結果。

由于SAMP 算法中步長是固定的,如果將步長s設置太大,則信號的真實稀疏性可能會被跳過,得到的結果可能并非全局最優(yōu)解,這將導致重建精度降低;如果步長s設置太小,迭代次數(shù)將會顯著增加,嚴重影響了算法的效率。因此,本文對SAMP 算法進行了兩方面的改進,即初始稀疏度K0估計和動態(tài)步長更新策略。

其中λ∈(0,1) 為步長變化率。本文實驗設定為η=0.1,λ=0.5。改進后的SAMP 算法完整步驟如下。

輸入傳感矩陣θ,相關度n,觀測向量{y1,y2,…,yn},迭代次數(shù)M,稀疏系數(shù)變化率閾值η,初始步長s,步長變化率λ,鄰接矩陣A。

輸出信號稀疏表示系數(shù)估計^x。

步驟1根據(jù)式(10)構造時空傳感矩陣A,y=(y1,y2,…,yn)T。

步驟2初始化:g=ATy,F0=?,K0=1。

步驟3取g中K0個最大值的索引組成F0。

步驟4如果,則K0=K0+1,重復步驟4;否則,到步驟5。

步驟5初始化:^x=0,r0=y,I=K0,k=1。

4 實驗結果和分析

4.1 量化指標

本文在真實地震數(shù)據(jù)進行量化實驗,來驗證本文算法的可行性和有效性。衡量數(shù)據(jù)重建效果的量化指標是信噪比(SNR)以及峰值信噪比(PSNR)。

其中,f是原始數(shù)據(jù),是重建的數(shù)據(jù),MSE是原始數(shù)據(jù)和重建數(shù)據(jù)的均方誤差。

4.2 稀疏變換矩陣對地震數(shù)據(jù)重建的影響

本文中的數(shù)據(jù)集來自實際中收集的完整地震數(shù)據(jù)。為了進行實驗分析和驗證,將其設為601 ×626切片數(shù)據(jù),用作原始數(shù)據(jù)集,如圖3(左)所示為原始地震切片數(shù)據(jù)。為了驗證重建算法的有效性,對其進行了50%高斯隨機抽樣,如圖3(右)所示。為了更好地顯示重建結果,計算重建前后的SNR 和PSNR,計算結果如表1 所示。

圖3 原始地震數(shù)據(jù)和50%采樣地震數(shù)據(jù)

表1 重建前后SNR 和PSNR 對比

為了驗證本文中使用的K-SVD 詞典學習在地震數(shù)據(jù)上的稀疏重建能力,設置了一種生成式壓縮感知(GCS)[10]算法進行比較。在數(shù)據(jù)恢復之前,使用K-SVD 算法訓練得到601 ×1052 的超完備字典,然后將學習到的超完備字典用作變換基,利用本文提出的迭代算法對50%采樣數(shù)據(jù)進行重建,重建結果如圖4(右)所示。使用生成式壓縮感知算法獲得的重建結果如圖4(左)所示。

圖4 使用GCS(左)或使用超完備字典(右)作為轉換基礎的重建結果

從圖5 的直觀重建結果可以看出,盡管生成式壓縮感知的重建結果可以獲得具有良好圖像質量的地震數(shù)據(jù),但是重建的地震數(shù)據(jù)與原始數(shù)據(jù)中的高頻信息并不十分吻合。這反映了基于神經(jīng)網(wǎng)絡的生成式壓縮感知重建方法的不確定性,不適用于地震數(shù)據(jù)重建。另一方面,不同地貌的地震數(shù)據(jù)差異很大,基于網(wǎng)絡訓練的模型無法普遍使用,導致該方法在地震數(shù)據(jù)重建領域的局限性。

為了驗證K-SVD 字典學習方法在地震數(shù)據(jù)上優(yōu)于傳統(tǒng)的稀疏變換方法,本文使用地震數(shù)據(jù)重建中常用的4 種常見稀疏變換矩陣進行數(shù)據(jù)重建作為對比實驗,它們分別是傅立葉變換、離散余弦變換、小波變換和曲波變換。表2 列出了每個稀疏基的重建結果的平均SNR 和PSNR。從表中可以看出,傅立葉變換矩陣的重建精度最差,而超完備字典是其中最佳的稀疏基。

表2 5 種稀疏基重建地震數(shù)據(jù)的SNR 和PSNR 對比

通過以上實驗證明,本文選擇的K-SVD 字典學習方法對地震數(shù)據(jù)的稀疏重建能力較強。

4.3 迭代算法對地震數(shù)據(jù)重建的影響

為了比較不同迭代算法之間的性能差異,使用超完備字典作為數(shù)據(jù)重建的稀疏變換基,并使用不同算法重建50%的采樣數(shù)據(jù)10 次并計算平均值。重建效果和運行時間如表3 所示。

表3 不同迭代算法對比

從實驗結果的分析可以看出,OMP 算法在重建精度及執(zhí)行效率方面均不佳。SP、CoSaMP 和IHT算法的運行時間較短,但重建精度相對較低,且需要設置稀疏性。IRLS 和SAMP 不需要設置稀疏度,并且具有更好的重建效果,但是運行時間相對較長。本文方法具有最佳的重建效果。并且與IRLS 和SAMP 算法相比,大大縮短了運行時間。

本文增設不同采樣率實驗,基于10%～70%樣本數(shù)據(jù)進行了重建。其中OMP、SP、CoSaMP、IHT 設置固定稀疏度K=50,SAMP 算法設置步長s=5。不同算法的重建效果,見圖5 和圖6。

圖5 不同采樣率下的SNR 對照圖

圖6 不同采樣率下的PSNR 對照圖

從實驗結果的分析來看,需要設置稀疏度的重建算法(OMP、IHT、SP、CoSaMP)無法很好地處理不同采樣率情況。而無需設置稀疏度的算法(IRLS,SAMP 和本文方法)在多種采樣率情況下效果相對較好,且高采樣率情況重建效果更佳。但是,SAMP算法需要設置步長,如果采樣率較低且步長設置得太高,則會導致跳過最準確的稀疏度。當采樣率較高時,如果步長較小,則將導致過多的迭代。在實驗中,本文方法和IRLS 可以隨著采樣率的增加而保持穩(wěn)定的增長,但是本文方法所需的時間比IRLS 所需的時間要短得多。

4.4 基于微電阻率成像數(shù)據(jù)的重建實驗

為了更好地驗證算法的泛化能力,增加了基于微電阻率成像數(shù)據(jù)的重建實驗。微電阻率成像數(shù)據(jù)也具有空間相關性的特征,但不同于地震數(shù)據(jù)的不規(guī)則特征。微電阻率成像數(shù)據(jù)的不完整部分相對較寬,具有一定的規(guī)律性。為了證明本文中的算法也可以解決此類問題,本文使用多種算法基于真實的微電阻率成像數(shù)據(jù)進行重建,原始數(shù)據(jù)大小為360×1000。重建結果顯示在圖7 中。對局部進行放大,顯示結果如圖8 所示。

圖7 原始電成像數(shù)據(jù)及不同算法重建結果

圖8 局部放大后的結果

從恢復效果可以看出,除去本文算法,其他算法的重建結果都有不同程度的“失真”。由于OMP、SP、CoSaMP 和IHT 算法需要預設稀疏度,稀疏度設置不佳,重建結果會稍微模糊。在實際生產過程中,具體的稀疏度通常未知,因此這些算法不適用于該類問題。IRLS 和SAMP 算法重建結果細節(jié)較好,但是所需的計算時間太長。而且,這6 種重建方法重建結果的缺陷部位均有明顯的重建痕跡。與其他算法相比,本文算法利用了空間相關信息,使重建結果細節(jié)更清晰,水平過渡更平滑,并且相比IRLS 和SAMP 大大減少了操作時間。證明了本文算法具有較強的稀疏重建能力和泛化能力。

5 結論

本文提出了一種基于時空約束壓縮感知的地震數(shù)據(jù)重構算法。該算法對傳統(tǒng)的壓縮感知理論模型進行了修改,增加了時空約束,并使用K-SVD 算法基于現(xiàn)有地震數(shù)據(jù)訓練并獲得超完備字典來代替?zhèn)鹘y(tǒng)變換。然后,對SAMP 算法進行了改進,提出初始稀疏估計和可變步長更新策略,從而大幅減少了迭代次數(shù),提高了重構精度。將生成式壓縮感知和傳統(tǒng)壓縮感知方法作為本文算法的對照實驗,并在真實地震數(shù)據(jù)及微電阻率成像數(shù)據(jù)上進行實驗,結果表明了本文提出的重建方法相比其他方法具有較強的重建能力。此外,本文方法與其他方法一樣,在低采樣率情況下無法很好地重建。生成式壓縮感知方法在本文實驗中沒有能夠展示出很好的結果,這是由于該方法約束不足而具有很強的不確定性,仍然將基于神經(jīng)網(wǎng)絡的重建方法作為未來的研究點,這可能將會是實現(xiàn)基于更低采樣率重建的新突破口。