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

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

中國石油大學(華東)地球科學與技術學院， 山東青島 266580

0 引言

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

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

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

本文提出了一種基于譜投影梯度組稀疏約束的全波形反演方法，利用波形信息對震源參數進行反演.針對震源空間信息稀疏而時間信息連續的特點，設計了組稀疏約束的目標函數，并采用譜投影梯度算法進行最優化反演.基于組稀疏約束的微地震震源參數譜投影梯度反演能夠獲得高精度的震源位置信息與子波信息.

1 原理

1.1 震源參數反演算法與梯度求取

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

(1)

其中u是壓力波場，v是介質速度，S(x,z,t)是時空項震源函數，可以拆分為S(x,z,t)=f(x,z)w(t)，其中f(x,z)表示純空間項，w(t)是純時間項或源子波.空間項可以對應不同空間位置的空間分布源或點源，本文的研究對象為點源.

全波形反演通過模型參數的更新來減小模擬數據和觀測數據之間的偏差，在更新過程中直接計算目標函數關于參數的導數(敏感核函數)時計算量較大，伴隨狀態法可以避免對該導數的直接計算，通過引入伴隨波場能夠計算目標函數關于模型參數的導數(Plessix，2006).傳統的目標函數使用L2范數定義：

(2)

其中uobs為觀測記錄，ucal為模擬記錄，使得接收點處的正演記錄和觀測波場之間的偏差最小化.模型向量m是時空源.

根據目標函數和聲波波動方程可以將伴隨狀態方程表示為：

(3)

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

目標函數關于震源參數的梯度可以表示為：

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

(4)

所以，震源的梯度就是伴隨波場，即由模擬波場和觀測波場的殘差進行波場反傳獲得.

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

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

震源參數反演求解組稀疏約束問題可以表示為：

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

(5)

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

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

(6)

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

(7)

在無約束的情況下，Friedlander等(1999)證明了在一定條件下，SPG對嚴格凸二次型具有超線性收斂速度.本文所采用的約束條件為L2,1范數的組稀疏約束，故可以將迭代投影到可行集合{x|‖x‖2,1<τ}：

(8)

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

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

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

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

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

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

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

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

2 震源參數反演測試分析

2.1 譜投影梯度反演與逆時定位對比

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

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

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

2.2 算法參數影響效果分析

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

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

圖2 不同τ值條件下目標函數在迭代過程中的下降曲線

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

圖3 采用不同參數值τ的震源位置和子波的反演結果

2.3 噪聲對反演結果的影響

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

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

圖4 譜投影梯度反演的抗噪性測試

2.4 子波頻率對反演結果的影響

在實際的水力壓裂過程中，微地震信號的傳播會受到地層介質參數以及注入流體參數等的影響，微地震信號通常是寬頻帶數據，震源頻率會影響定位的精度，下面對不同主頻頻率的震源子波進行測試，分析子波主頻對于反演結果的影響.

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

圖5 對不同主頻震源子波的反演結果

2.5 離散分布源反演結果

通過壓裂設備對地下巖層進行水力壓裂過程中會產生一系列微地震事件，微地震的震源之間通常距離很近，并且震源的激發時刻也不相同，同時在時間和空間上存在差異.為了解實際的壓裂效果，需要對一系列微震源進行監測.不同時刻激發的源信號之間存在串擾，通常會導致在計算結果中難以分辨鄰近震源，影響微地震監測結果.

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

圖6 不同時刻激發多震源的反演結果測試

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

3 背景模型影響效果分析

3.1 速度模型的不確定性對反演結果的影響

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

為了測試速度模型不同平滑程度對反演結果的影響，定義平滑度來表示數據平滑的程度，表示為：

(9)

其中，f表示處理前的信號，fs表示平滑后的信號，n表示數據所有元素的個數，i表示數據中的單個元素.

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

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

圖8 速度模型擾動對反演結果的影響

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

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

反演采用的速度模型需要比較接近實際地下介質參數，由于該方法只對震源時空參數進行反演，速度模型誤差會導致震源參數反演的誤差.

3.2 速度模型反演優化

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

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

圖10 利用反演后的速度模型獲得震源參數反演結果

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

4 結論

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

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

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

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