基于組稀疏約束的微地震震源參數(shù)譜投影梯度反演

唐杰， 劉英昌， 李聰， 高翔， 孫成禹

中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 山東青島 266580

0 引言

在過(guò)去的幾十年里，非常規(guī)油氣藏受到越來(lái)越多的關(guān)注.在以水力壓裂的方式對(duì)非常規(guī)油氣藏進(jìn)行增產(chǎn)過(guò)程中，儲(chǔ)層巖石破裂會(huì)發(fā)出微地震信號(hào)(Rutledge and Phillips，2003；Maxwell，2014).這些信號(hào)攜帶著豐富的震源參數(shù)信息，有助于分析儲(chǔ)層內(nèi)部的活動(dòng)斷裂網(wǎng)絡(luò)(李政等，2019).確定微地震事件時(shí)空演化規(guī)律已成為勘探開(kāi)發(fā)地震學(xué)的一個(gè)重要領(lǐng)域，地表布設(shè)檢波器由于操作簡(jiǎn)便且可選擇寬孔徑分布，因此被廣泛使用(Duncan and Eisner，2010).但是與沿監(jiān)測(cè)井記錄的數(shù)據(jù)相比，微地震信號(hào)在從地下傳播到地表的過(guò)程中能量會(huì)衰減(唐杰等，2021)，地表記錄的微地震數(shù)據(jù)的信噪比(Signal-to-Noise Ratio，SNR)較低(Maxwell and Urbancic，2001; Forghani-Arani et al.，2012)，增加了在含噪聲數(shù)據(jù)中拾取微地震信號(hào)的難度.地表微震資料的低信噪比給震源空間定位和發(fā)震時(shí)刻估計(jì)帶來(lái)了很大的挑戰(zhàn).

傳統(tǒng)基于P波和S波走時(shí)(Eisner et al.，2009；Zhebel and Eisner，2015)的走時(shí)反演等微地震事件定位方法存在對(duì)噪聲敏感的缺點(diǎn).Gajewski和Tessmer(2005)采用波動(dòng)方程逆時(shí)偏移(McMechan，1982)的思想對(duì)被動(dòng)源進(jìn)行逆時(shí)成像(Time-Reverse Imaging，TRI).Artman等(2010)使用自相關(guān)逆時(shí)成像條件生成震源位置的圖像，Nakata和Beroza(2016)和Tang等(2020)分別對(duì)逆時(shí)定位成像條件進(jìn)行了改進(jìn)，但是該方法受到速度模型不確定性的影響，存在聚焦不佳的問(wèn)題，并且地表排列的檢波器還會(huì)受到接收孔徑和采樣時(shí)間等因素的影響(Bazargani and Snieder，2016).

全波形反演(Full Waveform Inversion，F(xiàn)WI)方法是一種數(shù)據(jù)驅(qū)動(dòng)的模型參數(shù)估計(jì)方法.該方法利用記錄到的波形信息，構(gòu)建目標(biāo)函數(shù)并利用最優(yōu)化方法擬合合成波形和觀測(cè)波形來(lái)反演模型參數(shù)，其對(duì)初至拾取的依賴性較小，因此對(duì)噪聲有較強(qiáng)的魯棒性(Tarantola，1984；Plessix，2006).全波形反演主要用于速度模型構(gòu)建，微地震記錄中包含的信息也可以用于反演震源參數(shù)(Virieux and Operto，2009).Wu和McMechan(1996)在假設(shè)震源時(shí)間函數(shù)、介質(zhì)速度及密度等參數(shù)已知的條件下，通過(guò)FWI來(lái)確定震源空間位置、發(fā)震時(shí)刻以及與雙力偶相關(guān)的振幅和角度等參數(shù)，實(shí)現(xiàn)了雙力偶源參數(shù)的反演.Michel和Tsvankin(2014，2017)使用全波形反演獲得了垂直橫觀各向同性介質(zhì)中微地震震源的矩張量和空間位置等參數(shù)，這種算法要求對(duì)震源數(shù)量及其位置有著良好的初始估計(jì).Kaderli等(2015)采用了交替反演震源位置和子波參數(shù)的全波形反演方法，成功獲得了震源位置和震源子波，然而在兩個(gè)參數(shù)聯(lián)合反演過(guò)程中會(huì)出現(xiàn)局部極小值，影響反演效果.Sun等(2016)在聯(lián)合反演震源參數(shù)和介質(zhì)速度時(shí)提出將局部歸一化算子作為稀疏加權(quán)函數(shù)來(lái)約束模型空間，提高了反演收斂速度.Wang和Alkhalifah(2018)研究了基于互相關(guān)目標(biāo)函數(shù)的FWI方法以反演震源位置、震源函數(shù)和速度模型.Sharan等(2016，2017)提出基于稀疏度提升方法來(lái)反演微震震源波場(chǎng)，獲得震源的位置和震源時(shí)間函數(shù)，能夠較好地估計(jì)源時(shí)間函數(shù).Rodriguez(2018)使用了一種在反演中可以避免各參數(shù)梯度計(jì)算的粒子群算法實(shí)現(xiàn)對(duì)震源位置、發(fā)震時(shí)刻及矩張量的反演.Shekar和Sethi(2019)研究了一種稀疏約束算法，不需要對(duì)源的數(shù)量和類型進(jìn)行假設(shè)，并證明該算法具有很好的魯棒性.唐杰等(2022)提出了一種考慮地層各向異性衰減作用的微地震波場(chǎng)正演模擬與震源機(jī)制全波形反演方法，分析了震源區(qū)各向異性速度和衰減系數(shù)、噪聲及觀測(cè)系統(tǒng)對(duì)震源機(jī)制反演的影響.梯度投影法由Rosen(1960)提出，Goldfarb和Lapidus(1968)對(duì)其進(jìn)行了改進(jìn)，其基本思想為：當(dāng)?shù)c(diǎn)在可行域內(nèi)部時(shí)，取該點(diǎn)處的負(fù)梯度方向?yàn)榭尚邢陆捣较?；?dāng)?shù)c(diǎn)在可行域邊界上時(shí)，取該點(diǎn)處負(fù)梯度方向在可行域邊界上的投影產(chǎn)生一個(gè)可行下降方向.該方法的優(yōu)點(diǎn)是理論簡(jiǎn)單易操作，能有效解決大數(shù)據(jù)量問(wèn)題，但是在求取步長(zhǎng)時(shí)需要不斷投影來(lái)測(cè)試步長(zhǎng)，會(huì)產(chǎn)生巨大的計(jì)算量.Birgin等(2000)將非單調(diào)線性搜索和譜步長(zhǎng)的概念引入投影梯度法，提出了譜投影梯度法(Spectral Projection Gradient method，SPG)，其可以極大提高計(jì)算過(guò)程中的收斂速度.SPG已成功應(yīng)用于解決基追蹤、套索問(wèn)題以及二次回歸等問(wèn)題(van den Berg and Friedlander，2009；Dai and Fletcher，2005).

本文提出了一種基于譜投影梯度組稀疏約束的全波形反演方法，利用波形信息對(duì)震源參數(shù)進(jìn)行反演.針對(duì)震源空間信息稀疏而時(shí)間信息連續(xù)的特點(diǎn)，設(shè)計(jì)了組稀疏約束的目標(biāo)函數(shù)，并采用譜投影梯度算法進(jìn)行最優(yōu)化反演.基于組稀疏約束的微地震震源參數(shù)譜投影梯度反演能夠獲得高精度的震源位置信息與子波信息.

1 原理

1.1 震源參數(shù)反演算法與梯度求取

以二維聲波方法為例，震源波場(chǎng)隨著時(shí)間增加而向外傳播，在檢波器位置處記錄波場(chǎng)，震源波場(chǎng)可以采用如式(1)的聲波方程表示：

(1)

其中u是壓力波場(chǎng)，v是介質(zhì)速度，S(x,z,t)是時(shí)空項(xiàng)震源函數(shù)，可以拆分為S(x,z,t)=f(x,z)w(t)，其中f(x,z)表示純空間項(xiàng)，w(t)是純時(shí)間項(xiàng)或源子波.空間項(xiàng)可以對(duì)應(yīng)不同空間位置的空間分布源或點(diǎn)源，本文的研究對(duì)象為點(diǎn)源.

全波形反演通過(guò)模型參數(shù)的更新來(lái)減小模擬數(shù)據(jù)和觀測(cè)數(shù)據(jù)之間的偏差，在更新過(guò)程中直接計(jì)算目標(biāo)函數(shù)關(guān)于參數(shù)的導(dǎo)數(shù)(敏感核函數(shù))時(shí)計(jì)算量較大，伴隨狀態(tài)法可以避免對(duì)該導(dǎo)數(shù)的直接計(jì)算，通過(guò)引入伴隨波場(chǎng)能夠計(jì)算目標(biāo)函數(shù)關(guān)于模型參數(shù)的導(dǎo)數(shù)(Plessix，2006).傳統(tǒng)的目標(biāo)函數(shù)使用L2范數(shù)定義：

(2)

其中uobs為觀測(cè)記錄，ucal為模擬記錄，使得接收點(diǎn)處的正演記錄和觀測(cè)波場(chǎng)之間的偏差最小化.模型向量m是時(shí)空源.

根據(jù)目標(biāo)函數(shù)和聲波波動(dòng)方程可以將伴隨狀態(tài)方程表示為：

(3)

式中λ(x,z,t)表示伴隨波場(chǎng).

目標(biāo)函數(shù)關(guān)于震源參數(shù)的梯度可以表示為：

gs(x,z,t)=λ(x,z,t).

(4)

所以，震源的梯度就是伴隨波場(chǎng)，即由模擬波場(chǎng)和觀測(cè)波場(chǎng)的殘差進(jìn)行波場(chǎng)反傳獲得.

1.2 組稀疏約束譜投影梯度反演算法

流體注入引起的應(yīng)力變化導(dǎo)致裂縫的產(chǎn)生，這種破裂過(guò)程伴隨著裂縫尖端局部微地震能量的釋放，由于裂縫稀疏地散布在空間中，可以認(rèn)為微震事件的位置是稀疏的，但是微地震震源參數(shù)不僅包含了位置信息，還有時(shí)間信息，時(shí)間信息作為連續(xù)的信號(hào)并不具有稀疏性，所以對(duì)震源參數(shù)進(jìn)行反演時(shí)可以考慮進(jìn)行組稀疏約束(Shekar and Sethi，2019)，即在震源反演過(guò)程中對(duì)初始震源模型數(shù)組進(jìn)行分組，先對(duì)時(shí)間信息求L2范數(shù)，再在組水平進(jìn)行L1范數(shù)運(yùn)算，實(shí)現(xiàn)只對(duì)位置信息的組稀疏.

震源參數(shù)反演求解組稀疏約束問(wèn)題可以表示為：

‖m(X,t)‖2,1<τ,

(5)

本文選取譜投影梯度(SPG)算法(Birgin et al.，2000)對(duì)震源參數(shù)進(jìn)行反演.SPG使用非單調(diào)線搜索，在該算法中并不要求每一步都單調(diào)下降，只要求目標(biāo)函數(shù)值在有限步內(nèi)降低，這可能使目標(biāo)函數(shù)值出現(xiàn)暫時(shí)提高，但依然會(huì)保持整體下降的趨勢(shì).相較于單調(diào)線搜索法，非單調(diào)線搜索能夠避免進(jìn)入搜索狹區(qū)，進(jìn)而提高收斂速度.

對(duì)于給定的震源參數(shù)初始迭代點(diǎn)m0和正整數(shù)q，M，以及參數(shù)α>0， 0<ρ<1，非單調(diào)線搜索(Grippo et al.，1986)滿足：

(6)

SPG算法還使用了Barzilai和Borwein(1988)引入的譜步長(zhǎng)：

(7)

在無(wú)約束的情況下，F(xiàn)riedlander等(1999)證明了在一定條件下，SPG對(duì)嚴(yán)格凸二次型具有超線性收斂速度.本文所采用的約束條件為L(zhǎng)2,1范數(shù)的組稀疏約束，故可以將迭代投影到可行集合{x|‖x‖2,1<τ}：

(8)

其中τ為正則化參數(shù)，即范數(shù)球半徑，表示將震源參數(shù)m投影到L2,1組稀疏約束的范數(shù)球上.

在譜投影梯度算法中，投影操作是對(duì)目標(biāo)參數(shù)進(jìn)行反演的關(guān)鍵步驟，通過(guò)式(8)的組稀疏約束投影算子，在反演過(guò)程中可以依據(jù)震源參數(shù)時(shí)空特性對(duì)其進(jìn)行針對(duì)性投影，使位置的稀疏性和時(shí)間的連續(xù)性都得到較好的反演.

1.3 譜投影梯度反演算法流程

步驟2：設(shè)迭代次數(shù)為k，初始時(shí)k=0，給定初始參數(shù)m0，整數(shù)M>1,λ0=1,步長(zhǎng)0<αmin<αmax,α0∈[αmin,αmax]；

步驟3：判斷公式‖P(mk-gs(mk))-mk‖=0是否成立，成立則mk為最小點(diǎn)，可直接退出迭代；否則進(jìn)行步驟4；

步驟4：計(jì)算更新方向dk=P(mk-αkgs(mk))-mk； 步驟5：計(jì)算投影后的更新迭代點(diǎn)mp=mk+λkdk；

步驟7：mk+1=mk+λkdk，sk=mk+1-mk，yk=gs(mk+1)-gs(mk)，返回步驟3，如滿足迭代條件退出，輸出mk+1；否則進(jìn)行步驟8；

步驟9：mk+2=mk+1+αk+1gs(mk+1)，返回步驟3判斷是否滿足條件，如滿足輸出mk+2，否則繼續(xù)進(jìn)行迭代k=k+1.

2 震源參數(shù)反演測(cè)試分析

2.1 譜投影梯度反演與逆時(shí)定位對(duì)比

速度模型選取Marmousi模型的一部分，模型大小為1500 m×1250 m，震源在(750 m，625 m)位置處激發(fā)，震源函數(shù)為主頻為25 Hz的雷克子波，檢波器位于地表，間距30 m，共50個(gè)，震源及檢波器分布如圖1a所示.計(jì)算過(guò)程中選取的模型網(wǎng)格大小為5 m，記錄時(shí)間0.75 s，采用聲波方程有限差分正演模擬可以獲得微地震震源的波場(chǎng)記錄，如圖1b所示.

本文將微地震源視作點(diǎn)源，定位過(guò)程中拾取波場(chǎng)能量最大值處作為震源空間位置.圖1c為采用自相關(guān)成像條件的震源逆時(shí)定位結(jié)果，圖1d為采用譜投影梯度反演震源位置的結(jié)果，圖1e為采用譜投影梯度反演震源子波的結(jié)果.逆時(shí)定位方法和譜投影梯度法都可以獲取單個(gè)震源的空間位置，但是譜投影梯度反演出的震源位置成像結(jié)果更加精確，聚焦效果更好，聚焦效果會(huì)影響多個(gè)近鄰源的精確拾取.此外譜投影梯度反演算法能夠反演出震源子波，對(duì)比反演后的子波與真實(shí)子波，可以看出反演的子波與真實(shí)子波有較好的擬合，此外分析子波開(kāi)始振動(dòng)的時(shí)刻還可以獲得震源激發(fā)時(shí)刻.盡管采用譜投影梯度反演相較于逆時(shí)定位的計(jì)算時(shí)間有所增加，但該算法提升了計(jì)算精度，并且能夠獲得震源子波的時(shí)間信息.

圖1 (a) 速度模型，其中‘*’為震源，為檢波器； (b) 正演記錄； (c) 采用自相關(guān)成像條件的逆時(shí)定位結(jié)果； (d) 譜投影梯度反演的位置結(jié)果，虛線交點(diǎn)對(duì)應(yīng)真實(shí)震源位置； (e) 譜投影梯度反演的子波結(jié)果

2.2 算法參數(shù)影響效果分析

對(duì)于稀疏反演問(wèn)題，在采用基于L2,1范數(shù)球約束下的譜投影梯度反演算法時(shí)，其中一個(gè)關(guān)鍵參數(shù)是球半徑(τ)值，即正則化參數(shù)的選取，其在一定程度上會(huì)影響對(duì)于稀疏性的約束，對(duì)于算法中的參數(shù)可以根據(jù)測(cè)試分析進(jìn)行選取.

圖2為測(cè)試采用不同τ值條件下目標(biāo)函數(shù)在迭代過(guò)程中的下降曲線，目標(biāo)函數(shù)可以反映模擬結(jié)果與真實(shí)結(jié)果的擬合程度，目標(biāo)函數(shù)值越小，擬合程度越高.可以看出采用不同τ值的目標(biāo)函數(shù)在迭代20次后都能下降到較低的水平，當(dāng)τ值較小時(shí)，在迭代較少次數(shù)時(shí)反演就能獲得較好的擬合程度.

圖2 不同τ值條件下目標(biāo)函數(shù)在迭代過(guò)程中的下降曲線

圖3為采用不同τ值的震源位置與子波的反演結(jié)果.圖3a、3b和3c分別表示τ值為5、15、25時(shí)的震源位置的反演結(jié)果，可以看出不同τ值條件下都可以在真實(shí)震源處聚焦，但是當(dāng)τ較小時(shí)，反演得到結(jié)果在震源處聚焦的效果較好，這證明τ值越小時(shí)，對(duì)于反演過(guò)程中稀疏性約束的效果越好.圖3d為采用三種τ值的震源子波的反演結(jié)果與真實(shí)震源子波結(jié)果的對(duì)比圖，可以看出當(dāng)τ值為5時(shí)的震源子波反演擬合程度最高，隨著τ值的增加，雖然震源子波的整體波形變化不大，但是與真實(shí)子波相比，波形還是發(fā)生了細(xì)微的畸變.在實(shí)際資料處理過(guò)程追求計(jì)算速度將τ值設(shè)置在5左右，可以使目標(biāo)函數(shù)快速下降到較低水平.

圖3 采用不同參數(shù)值τ的震源位置和子波的反演結(jié)果

2.3 噪聲對(duì)反演結(jié)果的影響

實(shí)際采集的微地震記錄通常信噪比很低，所以反演定位中需要考慮噪聲對(duì)觀測(cè)記錄造成的影響，下面對(duì)低信噪比噪聲記錄進(jìn)行測(cè)試分析.

圖4為對(duì)譜投影梯度反演的抗噪性測(cè)試，圖4a為對(duì)圖1b的正演記錄中加入強(qiáng)隨機(jī)噪聲的微地震記錄，可以看出有效信號(hào)已經(jīng)淹沒(méi)在隨機(jī)噪聲中，無(wú)法有效地拾取初至到時(shí)信息，難以采用走時(shí)反演；圖4b為對(duì)圖4a的含噪微地震記錄進(jìn)行譜投影梯度反演的結(jié)果，可以看出反演的震源位置在真實(shí)震源附近聚焦，證明譜投影梯度法對(duì)震源位置的反演有較好的抗噪性；圖4c為對(duì)含噪數(shù)據(jù)進(jìn)行反演的震源子波結(jié)果，反演子波與真實(shí)子波的波形基本擬合.低信噪比資料使得反演子波出現(xiàn)了高頻振蕩，但其能量較弱，基本不影響后續(xù)子波的利用，說(shuō)明本文方法對(duì)隨機(jī)噪聲具有較好的魯棒性.

圖4 譜投影梯度反演的抗噪性測(cè)試

2.4 子波頻率對(duì)反演結(jié)果的影響

在實(shí)際的水力壓裂過(guò)程中，微地震信號(hào)的傳播會(huì)受到地層介質(zhì)參數(shù)以及注入流體參數(shù)等的影響，微地震信號(hào)通常是寬頻帶數(shù)據(jù)，震源頻率會(huì)影響定位的精度，下面對(duì)不同主頻頻率的震源子波進(jìn)行測(cè)試，分析子波主頻對(duì)于反演結(jié)果的影響.

圖5為選用不同頻率的震源子波進(jìn)行反演測(cè)試的結(jié)果.圖5a、5b、5c和5d為分別采用5 Hz、25 Hz、50 Hz、75 Hz主頻雷克子波的震源位置反演結(jié)果，可以看出當(dāng)主頻較低時(shí)(5 Hz)，定位結(jié)果出現(xiàn)了較大的偏差，隨著主頻頻率的增加，定位精度逐步提高，呈現(xiàn)出主頻越高，定位聚焦效果越集中的現(xiàn)象；圖5e為采用不同主頻的子波反演結(jié)果與真實(shí)子波的相似系數(shù)，通過(guò)相似系數(shù)可以評(píng)價(jià)反演子波與真實(shí)子波的擬合程度，可以看出相似系數(shù)在選用子波的主頻在5 Hz到25 Hz之間出現(xiàn)了顯著增加，并在25 Hz到75 Hz之間保持緩慢的增加，說(shuō)明在子波主頻頻率較低時(shí)反演的擬合程度較低，并隨著主頻的增加擬合程度呈現(xiàn)上升的趨勢(shì)，這主要是由于低頻信號(hào)波長(zhǎng)較長(zhǎng)，分辨率相對(duì)較低，而高頻信號(hào)的分辨率相對(duì)較高.由于微地震信號(hào)高頻會(huì)受到地層衰減以及噪聲等因素的影響，在微地震定位反演中可以采用分頻多尺度反演，先用低頻信息圈定大致震源范圍，再采用高頻信息精確定位.

圖5 對(duì)不同主頻震源子波的反演結(jié)果

2.5 離散分布源反演結(jié)果

通過(guò)壓裂設(shè)備對(duì)地下巖層進(jìn)行水力壓裂過(guò)程中會(huì)產(chǎn)生一系列微地震事件，微地震的震源之間通常距離很近，并且震源的激發(fā)時(shí)刻也不相同，同時(shí)在時(shí)間和空間上存在差異.為了解實(shí)際的壓裂效果，需要對(duì)一系列微震源進(jìn)行監(jiān)測(cè).不同時(shí)刻激發(fā)的源信號(hào)之間存在串?dāng)_，通常會(huì)導(dǎo)致在計(jì)算結(jié)果中難以分辨鄰近震源，影響微地震監(jiān)測(cè)結(jié)果.

圖6為對(duì)不同時(shí)刻激發(fā)的多震源進(jìn)行反演結(jié)果的測(cè)試.圖6a為真實(shí)速度模型以及震源位置，三個(gè)離散源的位置分別為(650 m, 500 m)、(700 m, 625 m)、(750 m, 700 m)；圖6b為采用譜投影梯度反演的震源位置結(jié)果，可以看出反演的震源位置與真實(shí)震源位置重合度較高，在真實(shí)震源處可以較好地聚焦，說(shuō)明該反演算法可以對(duì)多震源的位置進(jìn)行較好的反演；圖6c為針對(duì)不同時(shí)刻激發(fā)的微地震進(jìn)行反演獲得的震源子波反演結(jié)果與真實(shí)震源子波的對(duì)比圖，可以看出，本文方法對(duì)不同時(shí)刻激發(fā)的震源子波反演結(jié)果與真實(shí)震源有較高的擬合程度，近鄰震源帶來(lái)的串?dāng)_噪聲能量較低，說(shuō)明該方法能夠較好地反演不同時(shí)刻多震源的信息.

圖6 不同時(shí)刻激發(fā)多震源的反演結(jié)果測(cè)試

測(cè)試結(jié)果證明了本文所提出的反演算法的有效性，基于組稀疏約束的反演方法能夠獲得時(shí)空上滿足所期望特征的源函數(shù)，即空間稀疏性以及時(shí)間連續(xù)性.

3 背景模型影響效果分析

3.1 速度模型的不確定性對(duì)反演結(jié)果的影響

在采用譜投影梯度反演震源參數(shù)的過(guò)程中需要依賴背景介質(zhì)速度模型進(jìn)行計(jì)算，但是真實(shí)地下介質(zhì)模型通常存在很多的不確定性，速度偏大或者偏小都可能會(huì)對(duì)最后的反演結(jié)果產(chǎn)生影響，需要分析速度模型的不確定性對(duì)譜投影梯度反演結(jié)果的影響.

為了測(cè)試速度模型不同平滑程度對(duì)反演結(jié)果的影響，定義平滑度來(lái)表示數(shù)據(jù)平滑的程度，表示為：

(9)

其中，f表示處理前的信號(hào)，fs表示平滑后的信號(hào)，n表示數(shù)據(jù)所有元素的個(gè)數(shù)，i表示數(shù)據(jù)中的單個(gè)元素.

圖7a為平滑度為78的震源位置反演結(jié)果，可以看出反演位置與真實(shí)震源位置出現(xiàn)了偏差；圖7b為其對(duì)應(yīng)的子波反演結(jié)果，反演結(jié)果同樣存在偏差；圖7c為采用6種不同平滑程度的模型進(jìn)行譜投影梯度法反演震源位置結(jié)果與真實(shí)震源之間距離的絕對(duì)誤差值，圖7d為采用6種不同平滑度情況下反演的震源子波與真實(shí)子波之間的殘差，兩者的誤差都隨著平滑度的增加而上升，說(shuō)明平滑度會(huì)對(duì)震源子波的反演精度產(chǎn)生影響.圖8為速度擾動(dòng)對(duì)譜投影梯度反演結(jié)果的測(cè)試，圖8a為背景速度模型增大5%后的位置反演結(jié)果，可以看出速度增大后，定位較真實(shí)震源位置偏淺，主要是由于速度值偏大導(dǎo)致傳播時(shí)間減少；同理當(dāng)速度值減小時(shí)如圖8c中所顯示，定位結(jié)果偏深，符合速度模型擾動(dòng)所產(chǎn)生的誤差；圖8b為速度模型整體增大5%的子波反演結(jié)果與真實(shí)震源子波的對(duì)比，可以看出反演出的子波存在滯后的現(xiàn)象；圖8d為速度減小5%的子波反演結(jié)果與真實(shí)震源子波的對(duì)比圖，子波出現(xiàn)超前現(xiàn)象，說(shuō)明速度模型的正擾動(dòng)和負(fù)擾動(dòng)分別會(huì)對(duì)子波的反演造成滯后和超前現(xiàn)象.背景速度模型的準(zhǔn)確性對(duì)震源參數(shù)的反演具有重要作用.

圖7 譜投影梯度反演的速度不確定性測(cè)試

圖8 速度模型擾動(dòng)對(duì)反演結(jié)果的影響

與大多數(shù)FWI算法一樣，震源定位全波形反演算法對(duì)速度模型具有敏感性，要求速度模型接近真實(shí)速度模型.速度模型之間的巨大差異導(dǎo)致了周期跳變，即預(yù)測(cè)數(shù)據(jù)與觀測(cè)數(shù)據(jù)相差超過(guò)半個(gè)周期，可能會(huì)導(dǎo)致反演收斂到局部極小值.根據(jù)之前使用的真實(shí)速度模型和震源位置(圖1a)獲得的微地震數(shù)據(jù)，采用如圖9a所示的一維速度模型反演震源位置，由圖9b所示實(shí)際震源位置和真實(shí)震源位置存在差異，由于周期跳變的影響，真實(shí)震源子波和反演獲得的震源子波之間存在明顯差異(如圖9c所示).

圖9 采用一維初始速度模型反演結(jié)果

反演采用的速度模型需要比較接近實(shí)際地下介質(zhì)參數(shù)，由于該方法只對(duì)震源時(shí)空參數(shù)進(jìn)行反演，速度模型誤差會(huì)導(dǎo)致震源參數(shù)反演的誤差.

3.2 速度模型反演優(yōu)化

速度模型和微地震事件描述問(wèn)題緊密聯(lián)系在一起，微震震源位置的表征是在一個(gè)不完全已知的背景介質(zhì)模型下進(jìn)行的，震源參數(shù)反演的準(zhǔn)確度依賴于背景介質(zhì)參數(shù).在非常規(guī)油氣開(kāi)發(fā)過(guò)程中，儲(chǔ)層是在不斷改造、不斷發(fā)生變化的，導(dǎo)致微地震震源定位結(jié)果需依賴不斷發(fā)生變化的介質(zhì)參數(shù)模型(張曉林等，2013).用于被動(dòng)源反演的初始速度模型是從測(cè)井曲線或走時(shí)層析成像等方法得到的，大多數(shù)情況下介質(zhì)速度模型可能不準(zhǔn)確.FWI能夠提供更詳細(xì)的速度結(jié)構(gòu)，可以將壓裂射孔時(shí)產(chǎn)生的微地震數(shù)據(jù)或該壓裂區(qū)其他已知震源位置的能量較強(qiáng)的微地震數(shù)據(jù)用于FWI中更新速度模型.

由地表地震數(shù)據(jù)、測(cè)井信息等確定初始模型，然后利用后續(xù)產(chǎn)生的微震事件波形來(lái)更新背景速度模型，真實(shí)速度模型如圖1a所示，反演的初始模型采用一維速度模型如圖9a所示，采用震源區(qū)的6個(gè)已知的微地震記錄，震源分布如圖10a所示.采用編碼方式進(jìn)行速度模型反演(Krebs et al.，2009；Moghaddam et al.，2013)，迭代后獲得的反演速度模型如圖10b所示，由圖可知反演后的速度模型得到了一定程度的改善，成像效果與所采用的微地震源空間分布及其波路徑有關(guān)，通過(guò)速度模型反演能夠產(chǎn)生合理的速度更新.通過(guò)迭代后的速度模型反演獲得的震源定位結(jié)果如圖10c所示，與圖9b比較可知定位精度得到了明顯的提高.圖10d為反演獲得的震源子波，子波反演的精度也得到了明顯提高.圖10f為同樣采用震源區(qū)的6個(gè)已知的微地震記錄，采用震源編碼的方式對(duì)平滑初始速度模型(如圖10e所示)進(jìn)行速度反演更新獲得的速度模型，利用更新后的速度模型進(jìn)行震源參數(shù)反演獲得的震源定位結(jié)果如圖10g所示，震源子波如圖10h所示，對(duì)比圖7a和7b的平滑模型反演結(jié)果，利用更新后的速度模型反演獲得震源參數(shù)得到了進(jìn)一步的改善.

圖10 利用反演后的速度模型獲得震源參數(shù)反演結(jié)果

微地震震源數(shù)量、空間分布和觀測(cè)系統(tǒng)等都會(huì)影響速度模型的反演結(jié)果.微地震震源位于壓裂區(qū)，檢波器通常分布在地表或觀測(cè)井中，采集到的微地震數(shù)據(jù)中大部分能量是透射波.FWI中反射波通常傾向于更新模型的高波數(shù)部分，而透射波傾向于更新低波數(shù)部分(Alkhalifah，2016).微震震源之間距離相對(duì)較近會(huì)導(dǎo)致有限的照明.透射波反演特性以及有限的照明導(dǎo)致采用微地震數(shù)據(jù)反演速度模型的結(jié)果分辨率相對(duì)較低.盡管在模型反演測(cè)試時(shí)可以通過(guò)簡(jiǎn)單地在更廣的范圍內(nèi)設(shè)置更多的源來(lái)獲得較好的更新結(jié)果，但本文實(shí)驗(yàn)的震源設(shè)置相對(duì)更符合實(shí)際，反映了水力壓裂過(guò)程中實(shí)際發(fā)生的情況.雖然反演得到的速度模型不完善，但反演結(jié)果在運(yùn)動(dòng)學(xué)上是精確的，通過(guò)速度反演能夠產(chǎn)生合理的速度更新，進(jìn)而改善微地震震源參數(shù)的反演效果.此外還可以通過(guò)地表主動(dòng)源地震勘探數(shù)據(jù)結(jié)合微地震記錄進(jìn)一步改善速度模型參數(shù)，從而獲得精度更高的震源參數(shù)，本文不做進(jìn)一步展開(kāi).

4 結(jié)論

本文提出了一種基于譜投影梯度的組稀疏約束全波形反演微地震震源參數(shù)的方法，針對(duì)震源時(shí)間信息連續(xù)與空間位置稀疏的特點(diǎn)，在反演過(guò)程中采用L2,1范數(shù)進(jìn)行組稀疏約束，對(duì)反演正則化參數(shù)的選取進(jìn)行了測(cè)試，并對(duì)反演過(guò)程中的各種影響因素進(jìn)行了分析討論.研究結(jié)果表明：

(1)本文方法與逆時(shí)定位算法相比，空間源函數(shù)聚焦良好，與點(diǎn)源的真實(shí)位置非常接近.該方法可以同時(shí)獲得震源位置、發(fā)震時(shí)刻以及子波波形等信息，并且對(duì)不同時(shí)刻激發(fā)的多震源也具有較好的反演效果，基于組稀疏約束的反演方法能夠獲得時(shí)空上滿足所期望特征的源函數(shù).

(2)通過(guò)對(duì)反演過(guò)程中正則化參數(shù)進(jìn)行測(cè)試，得出正則化參數(shù)與源聚焦效果之間存在負(fù)相關(guān)，即源聚焦效果隨著正則化參數(shù)τ的減小而提升.該方法對(duì)噪聲具有魯棒性，高頻信號(hào)定位精度更高但易受噪聲影響，低頻信號(hào)受噪聲影響較小，但由于波長(zhǎng)較長(zhǎng)其定位分辨率相對(duì)較低，在微地震定位中可以采用分頻多尺度反演提高定位精度.

(3)微地震震源參數(shù)反演精度與背景速度模型的準(zhǔn)確性密切相關(guān)，微地震監(jiān)測(cè)中初始速度模型通常存在不確定性，影響了震源參數(shù)反演結(jié)果的可靠性.先通過(guò)壓裂數(shù)據(jù)進(jìn)行速度反演更新背景速度模型，然后利用更新后的速度參數(shù)進(jìn)行震源參數(shù)反演，能夠改善震源參數(shù)反演的結(jié)果，介質(zhì)與震源多參數(shù)聯(lián)合反演方法有待深入研究.