基于樣條插值與曲波變換壓縮感知的井下微地震監測數據重建

常 凱，張海江,b*,林 葉,b

(中國科學技術大學 a.地球與空間科學學院萬泰微地震實驗室，b.地震與地球內部物理實驗室，合肥 230026)

基于樣條插值與曲波變換壓縮感知的井下微地震監測數據重建

常 凱a，張海江a,b*,林 葉a,b

(中國科學技術大學 a.地球與空間科學學院萬泰微地震實驗室，b.地震與地球內部物理實驗室，合肥 230026)

對于非常規油氣開發水力壓裂井下微地震監測，由于井中布設的檢波器數量有限而導致空間采樣不夠，致使在進行偏移成像時產生空間假頻。為消除假頻現象，需要對空間欠采樣的數據進行內插。目前基于曲波變換(Curvelet)稀疏約束的地震道插值方法(即壓縮感知)已經逐漸被廣泛應用,但是對于井下微地震監測數據，由于檢波器數量過少，其在時間上的采樣點數要遠遠大于檢波器數。在這種情況下，曲波基函數在一定尺度下所能取到的方向信息是有限的，致使其基函數各向異性特征不能很好地發揮，因此插值效果不好。為解決該問題，我們提出了結合樣條插值的基于曲波變換稀疏約束地震道插值方法。對于井下微地震監測數據，首先利用樣條插值方法對每道數據的特定震相到達時間進行插值，在達到曲波基函數可以取到的方向后(如水平方向)，此時曲波基函數可以更稀疏地表示微地震數據體，然后再進行基于曲波變換稀疏約束的道間插值，并進行方向濾波，即可以達到較好的插值效果。我們基于一個井下微地震監測合成數據進行了測試，與數據加密前的地震偏移成像結果對比，內插之后數據偏移成像中的假頻得到了壓制。

井下微地震監測； 樣條插值； 曲波變換； 壓縮感知； 數據重建

0 引言

對于非常規油氣開發來說，例如頁巖氣和致密砂巖氣的開采，由于儲層滲透率低，必須采用水力壓裂技術使地層產生裂縫以提供油氣運移的通道。水力壓裂井下微地震監測，是一種有效的刻畫目標儲層裂縫分布的手段。井下微地震監測指的是把一串檢波器放在靠近壓裂井的監測井中，接收裂縫發育過程產生的微地震信號并對微地震事件進行定位或者對界面偏移成像。但是由于監測成本和技術等因素限制，井下微地震監測通常只布置十幾個檢波器[1]。在利用井下檢波器接收到的微地震波形進行偏移成像時，會由于道間距過大而出現空間假頻現象[2]。為解決此類問題，微地震數據道間插值成為抗空間假頻非常必要的一種手段[3]。

地震道集經過數學變換后可以用較少的基函數系數來表示[4]。Donoho等[5]提出了壓縮感知技術，在理論上解決了在數據不滿足Nyquist采樣定理的情況時，具有稀疏性質數據如何重建的問題。在此理論框架下，多種基函數被應用到地震數據恢復中，常用的有傅里葉基函數，小波基函數和曲波(Curvelet)基函數等[6-8]。其中曲波基函數在稀疏表示地震道集上表現出了相當好的各向異性，局部性和多尺度特性。

Herrmann等[9-10]在曲波框架下基于壓縮感知原理,對地震數據重建問題進行了深入的研究，取得了很好的效果。在國內，劉國昌等[11]將凸集投影(POCS)方法引入到基于曲波變換L1范數約束的地震道插值問題；張華等[12]采用了基于jilter采樣的曲波變換對三維地震數據進行了重建;馮飛等[13]將三維曲波變換和焦點變換結合起來，對三維地震數據進行插值和去噪，取得了良好的效果。然而，當已知地震道數很少，同時又有較長的時間采樣點數時(例如井下微地震監測數據)，若直接采用基于曲波變換稀疏約束的插值方法，由于曲波在具體尺度下的方向數是有限的，當同相軸不能符合其基函數的方向時，插值效果并不理想。

針對上述問題，筆者提出了一種樣條插值與基于曲波變換稀疏約束反演相結合的方法對井下微地震監測數據進行加密。這里首先分析了壓縮感知的基本原理并簡要說明曲波變換的方向性和多尺度特性；然后提出結合樣條插值和基于曲波變換稀疏約束數據重建的混合插值方法；最后基于一個井下微地震監測合成數據，比較了插值前后逆時偏移成像結果。

1 基于壓縮感知數據重建基本原理

地震數據插值可以看成由欠采樣數據到全部數據的一個數據重建反問題[14]。

d=M m

(1)

其中:d為觀測到的數據;M為采樣矩陣;m為想要恢復的全部地震數據,它可以通過數學變換W來進行稀疏表示，如式(2)所示：

x=Wm

(2)

通過對x的逆變換可以恢復m。如果數學變換W是一個正交變換，其逆變換與轉置等價：

m=WTx

(3)

因此，地震數據插值問題可以進一步轉化為稀疏約束優化問題，即得到的解一方面要擬合觀測數據，另外要保持其稀疏性，對應的優化問題的目標函數，如式(4)所示。

(4)

其中:G=MWT;λ為正則化參數;‖·‖0表示x中非零元素個數，即稀疏性。因為0范數約束方程求解時不可導，且為NP-hard問題，所以通常利用1范數替代[15]。最后得到最優化問題的目標函數為：

(5)

這里采用SPG_L1方法求解此泛函，其基本原理是將基追蹤(Basis Pursuit)問題轉化為一系列的最小絕對收縮和選擇(Least Absolute Shrinkage and Selection)問題進而求解[16]。因為曲波變換具有良好的多尺度特性和一定的方向特性，所以地震數據插值問題中的數學變換W可選為曲波變換。

在離散曲波變換中，由于其基函數的方向特性，空間不同方向的分量可以在曲波域中離散表示(如圖1)。對于二維函數f(n1,n2),0≤n1,n2≤n，其在笛卡爾坐標系下的曲波系數可以表示為：

(6)

圖1 數據在空間域與曲波域的表示Fig.1 Data representation in both space domain and curvelet domain(a)空間域數據；(b)對應于空間域數據的曲波系數

圖1(a)為空間域中兩個“同相軸”，在曲波域中可以表示為不同位置，不同尺度，不同角度稀疏的曲波系數(圖1(b))。由于一般情況下相鄰地震道“同相軸”不會在近垂直方向，因此可以通過對曲波系數按照方向特征濾波，濾掉在空間域中近垂直分量，達到提高插值效果的目的。

2 聯合樣條函數基于曲波變換的稀疏約束插值

基于壓縮感知技術插值的前提條件,是波形可以在變換域中稀疏表示。然而對于微震數據，由于井下檢波器數過少，致使曲波基函數無法直接稀疏表示全部微震數據。由于以上條件限制，我們將樣條插值與基于曲波變換的稀疏約束插值相結合。當樣條方向與曲波基函數所能稀疏表達的方向一致時，此時的曲波變換能夠使數據較好地稀疏表示，進而利用壓縮感知技術插值。其中三次樣條插值以其兼顧低次插值穩定性和高次插值光滑性的優點,已在工程中得到廣泛應用[18]。

對于空間欠采樣的井下微地震監測數據，為壓制微震成像時產生的空間假頻，采取以下技術流程來實現微震數據特定震相的重構。

1)對于高信噪比的數據，利用相鄰兩道的波形相似性并基于整體的背景速度，在微震波形中找到相同的震相，此震相視需要而定。

2)利用樣條插值擬合同相軸的形態，得到所有要插值的地震道上特定震相在時間軸上的位置。

3)基于樣條插值得到的相同震相的位置，對每個地震道在時間軸上截取相同的時間，然后使相同震相在曲波基函數能夠稀疏表達的方向上(水平方向)對齊。

4)對齊后的數據在曲波域中可以得到更稀疏的表示，滿足稀疏約束條件，進而利用壓縮感知技術進行插值并進行方向濾波。

5)將插值后的數據再返回到原來的位置上，最終實現微地震數據的特定震相的重構。

為了測試聯合樣條插值和曲波變換的數據重構效果，我們構建了一個井下微地震監測系統(圖2(a))。理論模型的背景P波速度為4 200 m/s，裂縫作為一個低速異常位于(1 000 m,2 200 m)附近，速度為3 570 m/s。34級檢波器(三角形)等間距地安裝在x=445 m的監測井處，深度范圍為975 m～1 965 m，空間間隔為30 m。在低速體和監測井之間(800 m,2 100 m)附近存在8個微地震事件(星形)。利用有限差分方法基于聲波方程計算了8個微地震事件的地震波場。圖2(b)顯示了一個微地震事件被34個檢波器接收到的地震合成波形。由圖2(b)可以看出，直達波以及對應于低速裂縫區域的散射震相。逆時偏移成像方法用于微地震成像時，如果是對微震位置定位，那么只需要將直達波反傳。因為后續散射波可認為由二次震源產生，在反傳時會在強散射點處加強，有可能在定位時造成一定程度的干擾；如果不考慮定位，比如已知位置的射孔事件，而是對裂縫等低速結構成像，那么就類似于傳統勘探手段中的逆時偏移方法。反傳時若含有直達波，同樣也會對最終成像結果產生干擾。對于實際數據來說(圖2(c))，我們可以較容易地拾取直達波震相，然而對于后續散射波的同相軸較難拾取，故而利用實際數據成像時選擇直達波后一段波形來進行加密。

圖3為該技術應用在模型上的算例。由于只對低速區界面成像，圖3(a)顯示了只包含后續散射震相的34個地震道。我們對原始數據進行均勻抽稀，得到了12個地震道，空間間距為90 m(圖3(b))，即稀疏度為66 %。首先拾取已有地震道相同震相最大振幅所對應的到時，并根據樣條插值計算丟失地震道的到時，然后確定數據重構的時間窗(圖3(c))，最后根據到時把時間窗向上移動進行展平(圖3(d))。在這個基礎上，對展平之后時間窗內的地震道根據公式(5)給出的曲波變換稀疏約束插值方法進行數據重構(圖3(e))，接著根據到時再把時間窗移回原來的位置(圖3(f))。為了比較，我們直接利用了基于曲波變換稀疏約束的內插方法進行數據重構(圖3(g))。為驗證該方法的普適性，用相同的技術路線測試了既包含直達波也包含后續散射波的欠采樣數據插值(圖3(h))，可以看出，直達波和散射波都具有較好的插值效果。

圖3 聯合樣條插值與曲波變換稀疏約束數據重構方法對一個微地震事件合成道集的測試Fig.3 Test on synthetic waveforms from a microseismic event for joint spline and curvelet-based sparsity constrained interpolation method(a)原始34道合成數據；(b)均勻抽稀后的12道原始數據； (c)利用樣條插值擬合散射震相的同相軸；(d)水平對齊后的數據；(e)對齊后應用曲波變換稀疏約束插值的數據；(f)最終插值結果；(g)直接利用曲波變換稀疏約束插值得到的數據；(h)包含了直達波和后續散射波的欠采樣波形用該技術插值效果

在僅有12道數據的情況下聯合樣條插值重構的地震數據，地震相位連續，效果良好(圖3(f))。與直接用曲波變換稀疏約束插值的結果相對比(圖3(g))，聯合方法的插值結果在相位連續性以及波形相似性更好。

為了進一步檢驗我們的方法對數據的恢復情況，計算了重構數據與正演合成數據的差值(圖4)。由圖4可以看出，聯合樣條插值的方法能夠更精確地恢復振幅(圖4(a)、圖4(c))。我們進一步選取了其中差值最大的一道數據進行波形對比(圖4(b)、圖4(d))。可以看出聯合了樣條插值的基于曲波變換稀疏約束的數據重構方法能夠很好地恢復相位、振幅。相比較而言，直接基于曲波變換的插值方法雖然也在一定程度上恢復了相位，但是存在一定程度的偏差，而且振幅恢復較差。

為定量說明重建后的數據和正演合成數據的差別，在這里把正演合成數據與重建數據的差值當作“噪聲”，把正演數據作為“信號”，引入信噪比的概念來定量評價插值數據的恢復質量[8]：

(7)

其中:m0為正演合成的34道數據；m1為重建數據；SNR單位為dB。經計算直接用曲波變換插值重構的數據信噪比僅為-1.8 dB，而聯合樣條插值和曲波變換的插值方法重構的數據信噪比為14.8 dB。

3 井下微地震合成數據加密前后成像效果對比

為檢驗該方法加密數據壓制空間假頻的有效性，將其應用于圖2所示井下微地震合成數據的偏移成像中。利用8個微地震事件的全部34道微震數據并基于聲波方程進行了逆時偏移成像(圖5(a))。其基本原理為從震源發出的正傳波場與井下檢波器所計算的逆推波場，在成像點處對傳播方向不同的波場應用互相關成像條件成像[19-20]。由圖5(a)可以看出，裂縫區域能夠很好地被刻畫出來。如果利用抽稀后的12道微震數據進行成像，雖然裂縫區域能夠較好地成像，但是在監測井附近存在較強的空間假頻現象(圖5(b))，這主要是由于微震波場空間采樣不夠所導致的。如果不聯合樣條插值，直接利用曲波變換進行數據重構得到34道微震數據，然后進行偏移成像(圖5(c))。由圖5(c)可以看出，雖然地震道進行了加密，但是由于插值效果較差，成像的結果依然顯示出與圖5(b)類似的空間假頻現象。與之不同的是，結合樣條插值基于曲波變換的數據重構方法能夠很好地恢復數據，而且成像結果中也不存在由于空間采樣不足而產生的空間假頻現象(圖5(d))。

通過對比圖5(a)～圖5(d)，我們還可以發現，加密后的數據與全部數據相比，在裂縫遠離監測井的界面成像效果稍有不足。加密后的數據和原始稀疏的數據相比，除了能夠更好地壓制假頻之外，對裂縫的成像效果并沒有太多形態上的差別。這說明基于樣條插值和曲波變換壓縮感知方法加密數據，并沒有帶來新的有效信息。

4 結論與討論

筆者提出了聯合樣條插值與曲波變換稀疏約束的數據重構方法，并用于壓制微地震偏移時由于空間采樣不足而產生的假頻?；诰挛⒌卣鸨O測合成數據測試，取得以下幾點認識：

1) 基于數學變換稀疏約束插值的前提，是地震道集圖像可以在變換域稀疏表示。對于曲波變換來說，只有當同相軸所表示的方向與其所對應的尺度下基函數所能刻畫的方向一致時，才可以稀疏表示全部地震道集，從而得到較好的插值效果。

2)聯合樣條插值與曲波變換稀疏約束的插值方法，可以針對在較少地震道的情況下，通過同相軸追蹤較好地進行地震數據重建。避免了當已知地震道較少時，曲波基函數對整個地震道難以進行稀疏表示的問題。

3)合成測試結果表明，針對空間采樣不足，用聯合方法加密微地震監測數據可以壓制偏移成像時出現的假頻現象，而并不帶來新的有效信息。

[1]張永華,陳祥,楊道慶,等.微地震監測技術在水平井壓裂中的應用[J].物探與化探,2013,37(6):1080-1084.ZHANG Y H,CHEN X，YANG D Q，et al.Application of microseismic monitoring technology in horizontal well fracturing [J].Geophysical and Geochemical Exploration,2013,37(6):1080-1084.(In Chinese)

[2]劉財,李鵬,劉洋,等.基于seislet變換的反假頻迭代數據插值方法[J].地球物理學報,2013 (5):1619-1627.LIU C,LI P,LIU Y,et al.Iterative data interpolation beyond aliasing using seislet transform.[J].Chinese Journal Geophysics,2013 (5):1619-1627.(In Chinese)

[3]林葉，張海江.水力壓裂裂縫微地震逆時偏移成像[C].2015年中國地球科學聯合學術年會——專題52:微地震監測與反演，2015:2273-2274.LIN Y,ZHANG H J.Reverse time migration imaging of microseismic data in hydraulic fracturing [C].Annual meeting of the Chinese Academy of Earth Science,2015.Special 52:microseismic monitoring and inversion，2015:2273-2274.(In Chinese)

[4]劉洋,王典,劉財.數學變換方法在地震勘探中的應用[J].吉林大學學報:地球科學版,2005,35(專輯) :1-8.LIU Y,WANG D,LIU C.The application of mathematical transform in seismic exploration [J].Journal of Jilin University:Earth Science Edition,2005,35 (Spec):1-8.(In Chinese)

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

[6]LUO T,LIU C,YANGX T,et al.Seismic data reconstruction based on iterative linear expansion of thresholds[J].Global Geology，2015,18(2) :127-133.

[7]崔興福,劉東奇,張關泉.小波變換實現地震道內插[J].石油地球物理勘探,2003,38(B11):93-97.CUI X F,LIU D Q,ZHANG G Q.Wavelet transform to achieve seismic interpolation [J].Oil Geophysical Prospecting,2003,38(B11):93-97.(In Chinese)

[8]馮飛,王德利,張亞紅,等.結合曲波變換的焦點變換在地震數據去噪和插值中的應用[J].物探與化探,2013,37(3):480-487.FENG F,WANG D L,ZHANG Y H,et al.The application of focal transform in combination with curvelet transform to seismic data denoising and interpolation[J].Geophysical and Geochemical Exploration,2013,37(3):480-487.(In Chinese)

[9]NAGHIZADEH M,SACCHI M D.Beyond alias hierarchical scale curvelet interpolation of regularly and irregularly sampled seismic data[J].Geophysics,2010,75(6):WB189-WB202.

[10]HERRMANN F J,HENNENFENT G.Non-parametric seismic data recovery with curveletframes[J].Geophysical Journal International,2008,173(1):233-248.

[11]劉國昌,陳小宏,郭志峰,等.基于Curvelet變換的缺失地震數據插值方法[J].石油地球物理勘探,2011,46(2):237-246.LIU G C,CHEN X H,GUO Z F,et al.Missing seismic data rebuilding by interpolation based on Curvelet transform [J].Oil Geophysical Prospecting,2011,46(2):237-246.(In Chinese)

[12]張華,陳小宏.基于 jitter 采樣和曲波變換的三維地震數據重建[J].地球物理學報,2013 (5):1637-1649.ZHANG H,CHEN X H.Seismic data reconstruction based on jittered sampling and curvelet transform [J].Chinese Journal Geophysics,2013 (5):1637-1649.(In Chinese)

[13]馮飛.結合稀疏變換的稀疏約束反演一次波估計研究[D].長春：吉林大學,2014.FENG F.A study of estimation of primaries by sparse inversion in combination with sparse transformation [D].Changchun:Jilin University,2014.(In Chinese)

[14]SHANG X F.Inverse scattering:theory and application to the imaging of the Earth's seismic discontinuities[D].Boston:Massachusetts Institute of Technology,2014.

[15]ASTER R C,BORCHERS B,THURBER C H.Parameter estimation and inverse problems[M].Waltham:Academic Press,2011.

[16]謝志鵬.壓縮感知的稀疏重構算法研究[D].南京：南京航空航天大學,2012.XIE Z P.Research on sparse reconstruction algorithms in compressive sensing [D].Nanjing:Nanjing University of Aeronautics and Astronautics,2012.(In Chinese)

[17]CANDES E,DEMANET L,DONOHO D,et al.Fast discrete curvelet transforms[J].Multiscale Modeling & Simulation,2006,5(3):861-899.

[18]周強,曹琳昱,曹中林,等.樣條曲線擬合初至波剩余靜校正方法[J].地球物理學進展,2015 (3):1329-1332.ZHOU Q,CAO L Y,CAO ZH L,et al.A method of residual static corrections based on fitting first break using spline function [J].Progress in Geophysics,2015(3):1329-1332.

[19]林葉.VSP 逆時偏移角道集域成像方法與應用[D].北京：中國地質大學 ,2014.LIN Y.VSP reverse time migration angle gather domain imaging method and application [D].Beijing:China University of geosciences,2014.(In Chinese)

[20]LIU F,ZHANG G,MORTON S,et al.An effective imaging condition for reverse-time migration using wavefield decomposition[J].Geophysics,2011,76:s29-s39.

Downhole microseismic monitoring data reconstruction based on spline interpolation and curvelet-based compressive sensing

CHANG Kaia，ZHANG Hai-jianga,b*，LIN Yea,b

(University Science and Technology of China a.Wantai Microseismic Lab of School of Earth and Space Sciences,b.Laboratory of Seismology and Physics of Earth's Interior,Hefei 230026,China)

For downhole microseismic monitoring of hydraulic fracturing of unconventional oil/gas reservoirs,the number of borehole receivers is limited and as a result the spatial sampling of seismic wavefield is sparse,which can lead to artifacts in the seismic migration image for downhole microseismic migration.To mitigate this spatial aliasing issue,spatial interpolation is needed.Currently,compressive sensing using the curvelet-based sparsity constraint has been gradually widely used for seismic trace interpolation.However,for the case of downhole microseismic monitoring,because the temporal sampling is much denser than the spatial sampling,the directional information obtained for curvelet bases at certain scales is limited,making it difficult to take full advantage of anisotropic features of curvelet bases.As a result,the result of seismic trace interpolation is poor.To solve this issue,we propose a hybrid method that combines spline interpolation and curvelet-based compressive sensing.First,spline interpolation method is used to calculate the reflection of interpolated traces using those of existing traces.Then the traces are shifted to a certain direction that could well be represented by the curvelet basis,such as the horizontal direction.After shifting the traces,the seismic gather is sparse in the curvelet domain and thus the curvelet-based spartisty constrained interpolation can then be used to more accurately reconstruct more traces in the space domain.We tested the hybrid method on a synthetic downhole microseismic dataset.It shows that the artifacts on seismic migration image using interpolated microseismic data are greatly reduced compare to the migration image from the original sparse data.

downhole microseismic monitoring; spline interpolation; curvelet transform; compressive sensing; data reconstruction

2015-11-13 改回日期：2015-12-24

國家自然科學基金項目(41274055)

常凱(1990-)，男，碩士，主要從事被動源地震成像，稀疏約束反演研究,E-mail：chao_13@mails.jlu.edu.cn。

*通信作者：張海江(1973-)，男，教授，主要從事不同尺度下地震成像算法研究和應用，E-mail：zhang11@ustc.edu.cn。

1001-1749(2016)06-0788-08

P 631.4