馬堅偉

(哈爾濱工業(yè)大學(xué)地球物理中心/數(shù)學(xué)系,黑龍江哈爾濱150001)

壓縮感知(Compressed Sensing,CS)描述的是可從高度不完備的線性測量中高精度重構(gòu)未知目標。創(chuàng)造性的把L1范數(shù)最小化和隨機矩陣理論有機結(jié)合,可得到稀疏信號重建的最佳效果。利用壓縮感知重構(gòu)一個連續(xù)信號,不再與香農(nóng)-奈奎斯特采樣定理所要求的“頻帶”有關(guān),而是與未知信號的“稀疏度”有關(guān)。CS由DONOHO,CANDS,TAO三位著名數(shù)學(xué)家提出[1-2],并入選了2007年美國十大科技進展,早期的主要幾篇文獻的引用已超過5×104次。

從數(shù)學(xué)上來理解,CS就是求解一個線性方程組Y=AX。Y是已知的觀察數(shù)據(jù),X是要重構(gòu)的目標。特別強調(diào)的是,這里重構(gòu)信號A是一個非方陣,其矩陣的行數(shù)遠少于列數(shù)。這相當于這個線性方程組的方程個數(shù)遠少于未知數(shù)的個數(shù)。顯然,這樣的方程組有無窮多個解,是一欠定反問題。但是CS告訴我們,如果先驗知道未知目標X有一定的稀疏性或者具有可壓縮性,且A的列之間具有一定的非相關(guān)性(比如隨機矩陣或勘探中的隨機接收器布點),那么可以用稀疏約束優(yōu)化的算法X=min‖Y-AX‖2+a‖TX‖0精確求解X。其中,T表示稀疏變換,零范數(shù)表示稀疏約束正則化。

滿足這個精確重構(gòu)所需要的測量數(shù)M可以表示為M=Slg(N/S),其中N表示待重構(gòu)信號的長度,S表示該信號的稀疏度。即便在有噪聲的情況下,即Y=AX+Noise,也可以進行高精度重構(gòu)。我們補充上面提到的兩個概念:稀疏和可壓縮性。稀疏是指一個信號如果其大部分振幅值為零,只有少數(shù)元素的振幅為非零,則稱該信號為稀疏信號;可壓縮性是指在某種變換域下的稀疏,即該信號在時空域下不稀疏,但可以找到一種稀疏變換,如傅里葉變換、小波變換、曲波變換等,使得該信號經(jīng)過變換后的系數(shù)只有少數(shù)非零值。我們可以從如下幾個角度來理解CS:從數(shù)學(xué)上講,CS本質(zhì)是降維,從低維空間去研究高維空間;從信號上講,CS本質(zhì)是采樣,從香農(nóng)采樣的頻率相關(guān)到稀疏度相關(guān);從工程上講,CS本質(zhì)是成本,從物理測量成本轉(zhuǎn)移到數(shù)學(xué)計算成本,用數(shù)學(xué)計算來彌補實際勘探的采集不足。

從20世紀40年代開始,統(tǒng)治信號和信息領(lǐng)域的一直是香農(nóng)-奈奎斯特采樣定理,即要想從離散信號中恢復(fù)連續(xù)信號,其離散信號的采樣率至少應(yīng)該是此連續(xù)信號最高頻率的兩倍?？梢哉f大部分相關(guān)信號采集的設(shè)備都是基于香農(nóng)定理設(shè)計的,如相機、雷達、核磁共振等醫(yī)學(xué)成像、通訊、信息等眾多工程領(lǐng)域的設(shè)備。CS理論出現(xiàn)后,其采樣率不再和信號的頻帶有關(guān),而是和信號的稀疏度有關(guān),從而打破了傳統(tǒng)理論的束縛。因此很多基于傳統(tǒng)理論設(shè)計的方法、軟件和設(shè)備都可以得到升級換代。2007年,萊斯大學(xué)提出的單像素相機就是一個很好的例子,不再去追求千萬像素的分辨率,而是基于壓縮感知理論用一個像素的時間序列成像就可重構(gòu)出高分辨率圖像。

壓縮感知是一個理論框架,不是單一的技術(shù),其三要素是隨機測量、目標的稀疏表示和稀疏促進的優(yōu)化重構(gòu)算法。地震勘探中,降低野外數(shù)據(jù)采集的成本而又能保證勘探精度是一個很重要的課題。反過來說,從現(xiàn)有采集的地震數(shù)據(jù)中能否挖掘出更多的寶貴信息來提高勘探精度？由于CS理論的出現(xiàn),這個問題由香農(nóng)采樣轉(zhuǎn)變到稀疏重構(gòu)的思路上,帶來了很多變革。CS的成功,在于改變了信號的采集方式。CS所涉及的技術(shù),在地震勘探的采集、處理、正演模擬、成像反演等方面都會帶來變革,但其核心問題之一還是地震數(shù)據(jù)和地質(zhì)目標的稀疏表示。

1 地震勘探中壓縮感知的最新研究進展

學(xué)術(shù)期刊Geophysics已多次組織了與壓縮感知地震勘探相關(guān)的專題,如2010年關(guān)于采樣和波場表示的專題、2015關(guān)于地震數(shù)據(jù)表示的專題。TheLeadingEdge也于2017年8月刊出了壓縮感知和高效采集的?？嚎s感知框架已經(jīng)滲入了傳統(tǒng)地震勘探的各個環(huán)節(jié),例如地震數(shù)據(jù)的采集、壓縮、重構(gòu)、模擬、偏移、反演等,并且在很多新的勘探技術(shù)中也發(fā)揮了作用,例如同時震源勘探、延時勘探等。不管應(yīng)用于何種技術(shù),地震數(shù)據(jù)的稀疏表示和隨機采樣始終是壓縮感知應(yīng)用的重要理論基礎(chǔ)。本章節(jié)首先介紹地震數(shù)據(jù)稀疏表示和隨機采樣的最新進展,然后分別從傳統(tǒng)勘探技術(shù)和新的勘探技術(shù)兩個方面介紹壓縮感知在地震勘探應(yīng)用中的最新進展。

1.1 稀疏變換

地震數(shù)據(jù)的稀疏表示按來源可分為3類。第一類是從圖像處理技術(shù)中借用已有的稀疏變換,如曲波變換[3-4]、剪切波變換[5-6]等等。曲波變換與剪切波變換均為多尺度多方向小波變換,適用于表示圖像中的線狀特征,因此也適用于表示地震數(shù)據(jù)的同相軸特征。第二類是利用字典學(xué)習(xí)方法從地震數(shù)據(jù)本身自適應(yīng)構(gòu)造稀疏變換,如K-均值奇異值分解(K-SVD)方法[7]、數(shù)據(jù)驅(qū)動緊框架方法(data-driven tight frame,DDT)[8]等等。這類方法通常在塊尺度下進行操作,適用于處理具有復(fù)雜結(jié)構(gòu)的地震數(shù)據(jù)。第三類是直接針對地震數(shù)據(jù)設(shè)計稀疏變換。物理小波變換[9]分別考慮了地震數(shù)據(jù)在時間上的小波狀和空間上的拋物狀特征,能夠?qū)ΟB前地震數(shù)據(jù)進行稀疏表示。地震小波變換[10-11]考慮了地震同相軸斜率變化的光滑性,實現(xiàn)了一種針對地震數(shù)據(jù)的提升小波。非對稱C-子波變換[12]考慮了地震子波的非對稱性,能夠在對地震數(shù)據(jù)稀疏表示的同時提取出相關(guān)信息,如走時、衰減參數(shù)、非對稱性參數(shù)等。

除了傳統(tǒng)意義上的向量的稀疏性,矩陣的低秩性(即奇異值的稀疏)也是一種廣義上的稀疏性,并且結(jié)合隨機采樣已經(jīng)應(yīng)用于地震數(shù)據(jù)重構(gòu)等[13]。此外,目前還出現(xiàn)了如拓展稀疏以及深度稀疏等新的稀疏性概念。拓展稀疏考慮了能量呈集中分布的情況;深度稀疏來源于深度學(xué)習(xí)技術(shù),即將稀疏基進一步進行稀疏表示,實現(xiàn)一個遞歸的稀疏表示過程。

1.2 隨機采樣

采樣過程的設(shè)計需要滿足一定的隨機性,同時滿足從高維到低維的采樣。隨機降維采樣在壓縮感知地震勘探如何體現(xiàn)呢？首先體現(xiàn)在觀測系統(tǒng)的設(shè)計上。我們可以采用不滿足奈奎斯特采樣定律的檢波器或者炮點分布來獲取采樣數(shù)據(jù),以此重構(gòu)完整數(shù)據(jù)。SHAHIDI等[14]引入了一種抖動(jitter)隨機采樣方法,即檢波器隨機分布在規(guī)則網(wǎng)格所限定的區(qū)域內(nèi)部,可控制相鄰檢波器的最大間距,避免數(shù)據(jù)大片丟失。NAGHIZADEH等[15]總結(jié)了幾種采樣方法,包括規(guī)則采樣、隨機采樣以及隨機加規(guī)則采樣,并研究了不同采樣的頻譜特征。NAGHIZADEH[16]還提出了一種雙編織三維地震數(shù)據(jù)采樣方法,這種采樣方法可以配合傅里葉方法進行數(shù)據(jù)重構(gòu)。HERRMANN[17]論述了利用隨機采樣甚至可從更少的數(shù)據(jù)中獲得更多的信息。國內(nèi)很多學(xué)者也做出了突出貢獻,例如陳生昌等[18]提出了幾種隨機與均勻結(jié)合的高效采集方法,并在數(shù)值試驗中取得較理想的效果。

除了觀測系統(tǒng)的下采樣設(shè)計,還可以利用隨機同時震源實現(xiàn)下采樣。具體來講,傳統(tǒng)勘探中不同炮是按次序激發(fā),現(xiàn)在可以同時激發(fā),減少了采樣時間,但采集中不同炮點的數(shù)據(jù)會混疊在一起。可以在同時激發(fā)中添加隨機性來解決這類問題,比如炮點激發(fā)時間有隨機延遲(炮編碼[19]),另外可以采用隨機的激發(fā)波形(波場模擬過程中可以實現(xiàn)[20])。這樣在共炮點域具有相關(guān)性的數(shù)據(jù),轉(zhuǎn)換到共檢波點域就是不相關(guān)的,進而可以進行分離。已有石油公司提出了有效的壓縮感知采集系統(tǒng),并已在實際勘探中進行了應(yīng)用測試[21-22]。

1.3 模擬、成像、混合震源、延時勘探

有限差分正演模擬的計算量取決于模型尺度以及觀測系統(tǒng)的分布,而不是取決于波場的復(fù)雜度。對于大尺度細密度的數(shù)據(jù)生成,即使波場非常簡單,仍需要大量的計算資源。HERRMANN等[23]提出了基于壓縮感知的同時全波形正演模擬,通過少量的獨立震源同時激發(fā),得到混疊波場,然后利用壓縮感知恢復(fù)算法得到分離后的波場。NEELAMANI等[20]考慮了波動方程格林函數(shù)的稀疏性,在震源同時激發(fā)時采用隨機帶限脈沖波形,利用壓縮感知理論和曲波變換來恢復(fù)原始波場。

VERA RODRIGUEZ等[24]實現(xiàn)了基于壓縮感知的微地震偏移成像技術(shù)。LI等[25]利用壓縮感知和隨機優(yōu)化技術(shù)實現(xiàn)了地震數(shù)據(jù)的快速全波形反演。作者使用了兩種采樣策略:隨機震源編碼和隨機震源位置。LI等[26]和李翔[27]利用稀疏約束的Gauss-Newton方法進行基于下采樣數(shù)據(jù)的全波形反演,減弱了下采樣所產(chǎn)生的人為噪聲,改善了反演結(jié)果。王華忠等[28]提出了基于特征波的成像方法,不是使用全波場而是使用壓縮感知框架提取出特征波場,保持走時特征,然后進行成像。

為了提高勘探效率,可以采用同時震源技術(shù)得到混疊的波場,然后利用壓縮感知恢復(fù)出沒有混疊的波場。ABMA等[29]利用傅里葉變換實現(xiàn)了混疊波場的分離,然后用于偏移成像,得到了與直接使用無混疊波場相似或更好的結(jié)果。原因是使用波場混疊技術(shù)可以在相同時間內(nèi)獲取更多炮點位置的波場記錄。LIN等[19]利用復(fù)雜的震源編碼技術(shù)以及曲波變換將不同的混合震源數(shù)據(jù)分離開來。

OGHENEKOHWO等[30]提出了基于壓縮感知的低成本延時地震勘探技術(shù)。不同時間的炮點布置在同一區(qū)域,但具體位置每次都不同,不同時間的勘探可以實現(xiàn)信息共享。另外,檢波器點是下采樣布置的,因此勘探成本較低。

2 討論與結(jié)論

CS如何在地震勘探開發(fā)中發(fā)揮作用？首先是針對數(shù)據(jù)采集,在CS理論框架下,可以用比傳統(tǒng)方法更少的采集成本實現(xiàn)同等(甚至更高)精度的反演。這個“更少”到底是多少,它與前面提到的反演目標的稀疏度有關(guān)。這個稀疏度(包括變換域的稀疏度),對一般的地震數(shù)據(jù)而言,是有一定的先驗知識;對實際數(shù)據(jù)的反演精度而言,只能有一個統(tǒng)計上的結(jié)論,即高概率重構(gòu)。此外,CS也可以深入挖掘現(xiàn)有勘探采集數(shù)據(jù)的信息,從而提高數(shù)據(jù)重構(gòu)(道加密)和偏移成像的反演精度。

如何才能用好CS？無非是在實際應(yīng)用中靈活把握CS的三要素:隨機采集(包括炮點和檢波器點兩方面的隨機)、目標的稀疏表達(明確具體問題的目標是高密數(shù)據(jù)還是地下介質(zhì)速度或其它物性參數(shù)？確定目標后,再找到或設(shè)計一個稀疏變換,進行目標的稀疏表達)和稀疏約束優(yōu)化重構(gòu)的快速算法。一般來說,重構(gòu)更高維的目標,需要用的采集數(shù)據(jù)(百分比)可更少。因為更高維的目標,有更多的先驗信息可以利用。CS在五維數(shù)據(jù)的插值重構(gòu)比在二維或三維數(shù)據(jù)的重構(gòu)更能發(fā)揮出優(yōu)勢,比如它有更多的斷層信息可利用,所以如何在高維目標的稀疏表達中更好利用結(jié)構(gòu)的幾何信息是一個很關(guān)鍵的課題。

經(jīng)過近十年的發(fā)展,壓縮感知理論已開始深入到地球物理勘探的方方面面,為降低勘探成本、提高采集效率、提高成像分辨率提供了新的理論支撐。相信在不久的將來,壓縮感知能夠真正應(yīng)用于實際油氣勘探,加速同時震源、頁巖油、深層勘探等非常規(guī)勘探發(fā)展。

壓縮感知可以理解為一種淺層的學(xué)習(xí),結(jié)合當前的深度學(xué)習(xí)技術(shù)是未來發(fā)展的一個方向。

致謝:感謝李幼銘、張捷、趙邦六、陳生昌、李成博等老師給予的巨大幫助,感謝課題組于四偉博士在此稿件的撰寫過程中做出的很多輔助工作。

[1] DONOHO D.Compressed sensing [J].IEEE Transactions on Information Theory,2006,52(4):1289-1306

[3] MA J,PLONKA G.The curvelet transform[J].IEEE Signal Processing Magazine,2010,27(2):118-133

[5] EASLEY G,LABATE D,LIM W Q.Sparse directional image representations using the discrete shearlet transform[J].Applied and Computational Harmonic Analysis,2008,25(1):25-46

[6] 馮飛,王征,劉成明,等.基于Shearlet變換稀疏約束地震數(shù)據(jù)重建[J].石油物探,2016,55(5):682-691

FENG F,WANG Z,LIU C M,et al.Seismic data reconstruction based on sparse constraint in the shearlet domain[J].Geophysical Prospecting for Petroleum,2016,55(5):682-691

[7] BECKOUCHE S,MA J.Simultaneous dictionary learning and denoising for seismic data[J].Geophysics,2014,79(3):A27-A31

[8] YU S,MA J,ZHANG X,et al.Interpolation and denoising of high-dimensional seismic data by learning a tight frame[J].Geophysics,2015,80(5):V119-V132

[9] ZHANG R,ULRYCH T J.Physical wavelet frame denoising [J].Geophysics,2003,68(1):225-231

[10] FOMEL S,LIU Y.Seislet transform and seislet frame [J].Geophysics,2010,75(3):V25-V38

[11] 劉財,李鵬,劉洋,等.基于seislet變換的反假頻迭代數(shù)據(jù)插值方法[J].地球物理學(xué)報,2013,56(5):1619-1627

LIU C,LI P,LIU Y,et al.Iterative data interpolation beyond aliasing using seislet transform[J].Chinese Journal of Geophysics,2013,56(5):1619-1627

[12] BOβMANN F,MA J.Asymmetric chirplet transform for sparse representation of seismic data[J].Geophysics,2015,80(6):WD89-WD100

[13] YANG Y,MA J,OSHER S.Seismic data reconstruction via matrix completion[J].Inverse Problems and Imaging,2013,7(4):1379-1392

[14] SHAHIDI R,TANG G,MA J,et al.Application of randomized sampling schemes to curvelet-based sparsity-promoting seismic data recovery [J].Geophysical Prospecting,2013,61(5):973-997

[15] NAGHIZADEH M,SACCHI M D.On sampling functions and Fourier reconstruction methods [J].Geophysics,2010,75(6):WB137-WB151

[16] NAGHIZADEH M.Double-weave 3D seismic acquisition,part 1:sampling and sparse Fourier reconstruction[J].Geophysics,2015,80(6):WD143-WD162

[17] HERRMANN F.Randomized sampling and sparsity:getting more information from fewer samples[J].Geophysics,2010,75(6):WB173-WB187

[18] 陳生昌,陳國新,王漢闖.稀疏性約束的地球物理數(shù)據(jù)高效采集方法初步研究[J].石油物探,2015,54(1):24-35

CHEN S C,CHEN G X,WANG H C.The preliminary study on high efficient acquisition of geophysical data with sparsity constraints [J].Geophysical Prospecting for Petroleum,2015,54(1):24-35

[19] LIN T T Y,HERRMANN F J.Designing simultaneous acquisitions with compressive sensing[J].Expanded Abstracts of 71stEAGE Annual Conference,2009:4-9

[20] NEELAMANI R N,KROHN C E,KREBS J R,et al.Efficient seismic forward modeling using simultaneous random sources and sparsity [J].Geophysics,2010,75(6):WB15-WB27

[21] LI C,MOSHER C,MORLEY L,et al.Joint source deblending and reconstruction for seismic data [J].Expanded Abstracts of 83rdAnnual Internat SEG Mtg,2013:82-87

[22] MOSHER C,LI C,JANISZEWSKI F,et al.Operational deployment of compressive sensing systems for seismic data acquisition [J].The Leading Edge,2017,36(8):661-669

[23] HERRMANN F J,ERLANGGA Y A,Lin T T.Compressive simultaneous full-waveform simulation[J].Geophysics,2009,74(4):A35-A40

[24] VERA RODRIGUEZ I,KAZEMI N.Compressive sensing imaging of microseismic events constrained by the sign-bit[J].Geophysics,2016,81(1):KS1-KS10

[25] LI X,ARAVKIN A Y,VAN LEEUWEN T,et al.Fast randomized full-waveform inversion with compressive sensing[J].Geophysics,2012,77(3):A13-A17

[26] LI X,ESSER E,HERRMANN F J.Modified Gauss-Newton full-waveform inversion explained—Why sparsity-promoting updates do matter [J].Geophysics,2016,81(3):R99-R112

[27] 李翔.基于壓縮感知技術(shù)的全波形反演[J].石油物探,2017,56(1):20-25

LI X.Full-waveform inversion from compressively recovered updates [J].Geophysical Prospecting for Petroleum,2017,56(1):20-25

[28] 王華忠,馮波,王雄文,等.特征波反演成像理論框架[J].石油物探,2017,56(1):38-49

WANG H Z,FENG B,WANG X W,et al.The theoretical framework of characteristic wave inversion imaging [J].Geophysical Prospecting for Petroleum,2017,56(1):38-49

[29] ABMA R,HOWE D,FOSTER M,et al.Independent simultaneous source acquisition and processing [J].Geophysics,2015,80(6):WD37-WD44

[30] OGHENEKOHWO F,WASON H,ESSER E,et al.Low-cost time-lapse seismic with distributed compressive sensing,part 1:exploiting common information among the vintages [J].Geophysics,2017,82(3):P1-P13