鄰近算法在微地震震源位置和一維速度模型同時反演中的應用

田宵 汪明軍 張雄 張偉 周立

1)江西省防震減災與工程地質災害探測工程研究中心(東華理工大學)，南昌 330013 2)南方科技大學，地球與空間科學系，廣東深圳 518055 3)中國地震局地震研究所，地震預警湖北省重點實驗室，武漢 430071

0 引言

近年來，隨著國內外頁巖氣開采的不斷深入，水力壓裂技術被廣泛地應用于非常規(guī)油氣藏的開發(fā)(張山等，2002；Warpinski，2009；Maxwell，2010)。水力壓裂過程中產生的微地震事件可用來監(jiān)測壓裂區(qū)域裂縫的發(fā)育情況和幾何形狀、評價壓裂效果以及監(jiān)測誘發(fā)地震，因此獲得準確的震源位置是微震監(jiān)測的基礎(Grechka，2010；Zimmer，2011)。微地震監(jiān)測主要包括地面監(jiān)測和井中監(jiān)測兩種方式(Eisner et al，2010；余洋洋等，2017，朱亞東洋等，2017)?？紤]到經濟成本和信號質量等因素，單井監(jiān)測是常見的微地震監(jiān)測方式，其主要優(yōu)點在于監(jiān)測到的微震信號信噪比相比地面監(jiān)測要高，且定位結果不受近地表的影響。

微地震監(jiān)測和常規(guī)的地震勘探不同，震源的位置和發(fā)震時刻均是未知的。因此，微地震定位方法通常借鑒于天然地震領域。目前常用的定位方法有很多種，如基于概率統(tǒng)計的貝葉斯方法(宋維琪等，2013；Zhang et al，2017)、基于走時的震源反演方法(Jones et al，2014；Zhou et al，2015；Tian et al，2016)以及基于波形的定位方法(王晨龍等，2013；Zheng et al，2016；Li et al，2018；田宵等，2020)。井下微震定位通常使用基于走時的定位方法，首先需要拾取微震信號的P波和S波初至，通過找到與初至曲線最吻合的空間位置來進行定位。該類方法主要包括利用絕對到時信息的Geiger經典定位方法(Geiger，1912)、利用走時差反演震源位置的相對定位算法(Waldhauser et al，2000；Guo et al，2017)以及全局的網格搜索方法(Jones et al，2014)。

一般來說，影響微地震單井定位結果的因素主要包括P波和S波到時、檢波器方位角和速度模型。其中，速度模型的準確性對定位結果至關重要。在井下監(jiān)測時，壓裂井和監(jiān)測井之間的距離為200～800m，通常將速度結構簡化為一維層狀模型。初始速度模型可從聲波測井曲線中提取，并利用已知震源位置的射孔事件進行速度校正(Warpinski et al，2005；尹陳等，2013；Jiang et al，2016)，然后進行震源定位。當沒有射孔事件時，速度模型可作為未知參數和震源位置同時反演(Li et al，2013；Tian et al，2017)。由于震源位置和速度模型之間有很強的耦合性，基于線性迭代的同時反演方法(本文簡稱為線性同時反演方法)通常依賴初始模型的選取，反演結果易陷入局部極小值。

近年來，全局反演算法因其不依賴于模型初始值的選擇而被廣泛應用于微地震速度模型反演。如李志偉等(2006)利用差異算法反演一維速度結構和震源參數；Jansky等(2010)采用鄰近算法和網格搜索方法分別更新一維速度模型和震源位置，并使用理論模型加以驗證；Tan等(2018)采用改進的鄰近算法和主臺站方法分別反演一維速度模型和震源位置。

鄰近算法(Neighbourhood Algorithm，簡稱NA)由Sambridge(1999)提出，其是一種解決非線性問題的直接搜索方法，用于在多維參數空間中尋找能夠擬合觀測數據的模型。該方法在非線性反演問題中應用較為廣泛，特別是觀測數據和未知參數之間的關系比較復雜的情況(Marson-Pidgeon et al，2000；Sambridge et al，2001)。為了研究鄰近算法在微地震單井監(jiān)測中的可行性，本文采用鄰近算法和網格搜索方法分別更新一維速度模型和震源位置(本文簡稱為鄰近算法-網格搜索)。首先通過設計不同層數的理論模型，研究速度模型的層數對鄰近算法精度的影響；然后將其應用于單井監(jiān)測的實際數據中，結合理論分析確定速度模型層數，并將速度模型和震源位置的反演結果與線性同時反演方法進行對比分析。

1 基本原理

1.1 利用網格搜索方法搜索震源位置

網格搜索方法是常用的微地震定位方法，具有簡單、快速等優(yōu)點。首先，需要對可能產生微地震事件的空間進行網格化，網格大小決定著搜索結果的精度以及搜索效率。然后利用射線追蹤方法計算每個網格點到檢波器的旅行時，最后用計算得到的旅行時與觀測到的走時數據進行匹配，誤差最小的網格點即最終確定的事件位置。當有大量地震事件需要進行定位時，為了節(jié)省計算量，對于給定的速度模型，可將計算的所有網格點到檢波器的旅行時保存成走時表，以提高運算效率。

需要搜索的震源參數包括震源位置和發(fā)震時刻，為了進一步提高搜索效率，本文使用P波和S波的走時差以及相鄰檢波器之間走時差消除發(fā)震時刻(Zhou et al，2016)，目標函數φ如下

(1)

其中，ns為震源的個數，mr為檢波器的個數，上標i和下標k、k′分別代表震源和檢波器，向量d為觀測的走時數據，G(m)為由模型m正演得到的計算走時。等式右邊第一項為每個微震事件的S波和P波的走時差，第二項為檢波器之間的P波走時差，第三項為檢波器之間的S波走時差。

1.2 利用鄰近算法搜索速度模型

鄰近算法是一種非導數搜索方法，該方法與遺傳算法、模擬退火屬于一類，其區(qū)別在于模擬退火和遺傳算法一般用來解決全局最優(yōu)解問題，而鄰近算法的目的是尋找參數空間中對數據擬合較好的區(qū)域進行優(yōu)先采樣的模型集合，而不是尋找單一的最優(yōu)模型，因此該方法可以更大程度地避免局部最優(yōu)解。

本文使用Sambridge開源的鄰近算法代碼(1)http://www.iearth.org.au/codes/NA/，在解決不同的實際問題時，一般需要調節(jié)2個參數(nm和nr)來滿足不同的采樣風格。nm為每次迭代時產生的模型數量，nr為每次迭代時重采樣的Voronoi單元的個數。根據Sambridge(1999)，一般nm的大小是參數維數的2倍，nr的范圍在2～nm/2之間，且要有足夠多的迭代次數。根據經驗，迭代次數可設置為nm的10～100倍。毋庸置疑，參數nm和nr的數值越大，就需要更多的迭代次數來尋找準確的解(尤其在搜索參數的數量較多時)。

聲波測井曲線可以較好地描繪出巖性分界面，層狀模型的厚度可直接從聲波測井曲線中提取。為了減少搜索參數的數量，不搜索速度模型的層厚。首先，使用鄰近算法產生若干個一維速度模型；然后對于每個速度模型，利用網格搜索方法進行震源定位并計算每個速度模型對應的走時殘差；再根據走時殘差產生新的速度模型；最后選擇合適的迭代次數、nm和nr來保證殘差函數收斂到最小值。詳細的算法流程如圖1 所示。

圖 1 鄰近算法-網格搜索反演速度模型和震源位置的算法流程

2 理論模型計算

一個單井監(jiān)測的二維理論模型如圖2(a)所示，包含4個微震事件(黑色圓圈)和24個檢波器(黑色三角形)。速度模型為8層的一維速度模型，鄰近算法的搜索范圍如圖2(b)中陰影區(qū)域所示。首先假設每層的縱橫波速度比相同，則搜索參數為6，包括5個P波速度和1個縱橫波速度比。

圖 2 理論測試模型示意圖(a)理論觀測系統(tǒng)，紅色框為網格搜索范圍；(b)一維水平層狀模型，灰色陰影區(qū)域為鄰近算法搜索范圍

采用Sambridge提供的鄰近算法代碼，選取4000個初始樣本，nm和nr分別為200和100，迭代次數為500，最后產生了104000個速度模型。網格搜索方法的橫向和縱向搜索步長均為2m，搜索范圍如圖2(a)中紅色框所示。反演結果表明在沒有噪音的情況下，全局反演方法可以得到真實的震源位置和速度模型。圖3(a)為鄰近算法-網格搜索方法搜索的前20000個速度模型對應的均方根(Root Mean Square，簡稱RMS)走時殘差，結果表明鄰近算法-網格搜索方法經過4000個初始采樣點后很快找到正確的采樣區(qū)域，走時殘差降至0.005s以下。最后僅需221次迭代即可搜索到正確的速度模型。

為了模擬實際數據的走時拾取誤差，對計算的旅行時添加1ms的隨機噪音。反演結果表明當數據有1ms噪音時，震源位置的平均定位誤差約16m，搜索過程的走時殘差如圖3(b)所示。通過對比發(fā)現(xiàn)，數據噪聲對RMS走時殘差有一定的影響，對搜索過程的收斂速度影響較小。為了更加清楚地展示搜索過程中產生的速度模型，圖4 給出了數據有噪音時鄰近算法-網格搜索方法搜索的前20000個速度模型，速度模型的顏色越接近紅色，代表RMS走時殘差越小。由圖4 可見，數據有誤差時，搜索的速度模型也有一定偏差。

圖 3 鄰近算法-網格搜索方法在沒有噪音(a)和有噪音(b)時的走時殘差分布

圖 4 數據有噪音時鄰近算法-網格搜索方法反演的P波(a)和S波(b)速度模型黑色虛線為真實速度值；藍色虛線為最終反演結果

以上理論測試中固定了縱橫波速度比，降低了搜索維度。但是實際情況的縱橫波速度比與泊松比有關，不同的巖石類型其縱橫波速比有所不同，除了一些特殊的巖石，縱橫波速度比的范圍在1.5～2.2之間。在處理實際數據時，很難將每層的縱橫波速度比當成固定值。對于一個n層的速度模型，搜索參數的個數應為2n。因此，討論搜索參數的數量對反演結果準確度和收斂性的影響是非常有必要的。圖5 給出了搜索參數的數量和RMS走時殘差的關系，為了更好地研究參數數量對反演精度的影響，圖中未對走時信息添加噪音。每個搜索參數均根據Sambridge(1999)給出的經驗調節(jié)nm、nr以及迭代次數。圖5 表明反演結果的RMS走時殘差隨著搜索參數的增加而急劇上漲。當搜索參數等于14時，對應的走時殘差量級為 10-4s。 因此，在利用鄰近算法-網格搜索方法反演一維速度模型時，應盡可能地減少速度模型的層數。

圖 5 不同的搜索參數數量對應的走時殘差

3 單井監(jiān)測的實際數據

將鄰近算法-網格搜索方法應用于單井監(jiān)測的水力壓裂實際數據中，該套數據為水平井壓裂、垂直井監(jiān)測。監(jiān)測井中包含了12個三分量檢波器，深度2800～3020m，間距20m。為了驗證鄰近算法-網格搜索方法的結果，對該套數據進行了線性同時反演，線性同時反演算法需要初始速度模型和初始震源位置。圖6(a)中灰色實線為聲波測井曲線，為了減少鄰近算法-網格搜索方法搜索參數的個數，本文從中提取了一個包括P波和S波的五層初始速度模型(黑色實線)，圖6(a)中黑色三角形為檢波器的深度位置?？紤]到檢波器和壓裂井所在區(qū)域，圖6(b)中灰色區(qū)域為鄰近算法-網格搜索方法搜索的速度范圍，搜索參數的個數為8。

圖 6 實際數據速度模型(a)測井曲線以及提取的初始速度模型；(b)鄰近算法-網格搜索方法搜索的速度范圍(灰色區(qū)域)

在第18個壓裂段選取40個信噪比較高的微震事件，并手動拾取了P波和S波到時?？紤]到一維速度模型和垂直監(jiān)測井的對稱性，將所有的地震事件投影到一個垂直平面(XZ)，將三維的定位問題簡化為二維，垂直平面包括一個垂直的Z坐標和一個水平的X坐標(事件到監(jiān)測井的距離)。隨后在二維平面上進行定位，并利用事件的方位角將定位結果反投影回三維空間，事件的方位角可由P波偏振分析確定。

分別用線性同時反演算法和鄰近算法-網格搜索反演震源位置和一維速度模型。對于該套數據，根據Sambridge(1999)的算法，選取1000個初始樣本，nm和nr分別為20和10，迭代次數為200，鄰近算法最終會產生5000個速度模型。圖7(a)為鄰近算法-網格搜索方法得到的所有速度模型對應的RMS走時殘差，結果表明200次迭代足以得到穩(wěn)定的結果，最優(yōu)的速度模型為第4947個，如圖7(b)中紅色虛線所示，對應的走時殘差為1.64×10-3s。表1 對比了不同的鄰近算法參數得到的RMS走時殘差，結果表明增加搜索樣本數或者更改速度模型的搜索范圍對反演結果的影響不大。

圖 7 鄰近算法-網格搜索反演結果(a)鄰近算法-網格搜索方法得到的速度模型對應的RMS走時殘差；(b)速度反演結果，黑色實線代表初始速度模型，紅色虛線為鄰近算法-網格搜索方法反演的結果，藍色虛線為以交替反演結果為初始值的線性同時反演結果

表 1不同的鄰近算法參數對應的RMS

對于線性同時反演方法，一般采用測井曲線提取的速度模型為初始速度模型(圖6(a)中黑色實線)，網格搜索結果為初始震源位置。采用不同方法對40個微震事件進行定位，結果如圖8 所示，其中包括X-Y平面和X-Z平面。進行水力壓裂作業(yè)時，微地震事件一般發(fā)生在壓裂段的附近，而網格搜索的定位結果(圖8 中綠色點)在深度上與壓裂井相差100m左右。以網格搜索的結果為初始值進行線性同時反演，得到的震源位置如圖8 中紫色點所示，對應的走時殘差為2.91×10-3s。從X-Z平面上來看，綠色點和紫色點的分布相似，均在壓裂井的下方。圖8 中紅色點為采用鄰近算法-網格搜索得到的定位結果，可以看出微震事件分布在壓裂井的兩側，反演結果較為合理。進一步對比鄰近算法-網格搜索方法(紅色點)和線性同時反演方法(綠色點)的定位結果和走時殘差，結果表明2種方法的定位結果在水平方向上相似，深度上差別較大。而線性同時反演的深度與初始位置相似，由此可推斷出，以網格搜索結果為初始值的線性同時反演方法可能陷入局部最優(yōu)解。此外，本文還采用Tian等(2018)的反演流程，首先利用交替迭代方法得到一個較好的初始解(包括震源位置和速度模型)，再采用線性同時反演方法，定位結果如圖8 中藍色點所示，速度反演結果為圖7(b)中藍色虛線，走時殘差為1.74×10-3s。對比鄰近算法-網格搜索方法(紅色點)和以交替迭代為初始值的線性同時反演(藍色點)結果，可以看出2種方法得到的定位結果均較為合理。線性同時反演方法容易陷入局部最優(yōu)解的主要原因為，微地震實際數據資料的信噪比較弱以及單井監(jiān)測的檢波器分布范圍較小，導致震源位置和速度模型之間存在較強的耦合性，求解空間存在多極值問題。因此，線性同時反演依賴初始值的選取，一個好的初始值有助于線性同時反演方法跳出局部最優(yōu)解。

圖 8 40個微震事件的定位結果綠色點為網格搜索的結果；紅色點為鄰近算法-網格搜索算法的結果；紫色點和藍色點分別為線性同時反演方法以網格搜索和交替迭代為初始值的結果

為了進一步評價速度模型和震源位置反演的精度，給出了2種方法反演的40個微震事件的走時擬合殘差分布直方圖，見圖9。走時擬合殘差為反演結果(震源位置和速度模型)正演的理論走時與拾取走時之間的偏差。圖9 中“同時反演1”為線性同時反演方法使用網格搜索的震源位置和測井速度為初始值時得到的走時殘差分布，而“同時反演2”為線性同時反演方法使用交替迭代結果為初始值時得到的走時殘差分布。通過對比，“同時反演1”、“同時反演2”和鄰近算法-網格搜索方法分別有49.8%、75.0%和74.8%的P波走時殘差分布在0～1ms范圍內，而3種方法的S波走時殘差分別有28.6%、56.1%和58.0%的事件分布在0～1ms范圍內。對比走時殘差的平均值，“同時反演1”、“同時反演2”和鄰近算法-網格搜索方法的P波走時殘差平均值分別為0.0012、7.65×10-4和7.62×10-4，而S波走時殘差的平均值分別為0.0021、0.0011和0.0011。因此，當線性同時反演使用較好的初始值時，得到的走時殘差分布與鄰近算法-網格搜索的結果大致相同。

圖 9 線性同時反演方法和鄰近算法-網格搜索方法的走時殘差直方圖(a)和(d)分別為線性同時反演方法使用網格搜索為初始值的P波和S波走時殘差分布；(b)和(e)分別為線性同時反演方法使用交替迭代結果為初始值的P波和S波走時殘差分布；(c)和(f)分別為鄰近算法-網格搜索方法的P波和S波走時殘差分布

上述結果表明，與線性同時反演方法相比，鄰近算法-網格搜索方法受初始模型影響較小，能夠最大限度地避免局部最優(yōu)解。雖然該方法能夠從局部最優(yōu)值中跳出，但其缺點為計算量較大。對比線性同時反演和鄰近算法-網格搜索方法在不同微地震規(guī)模下的計算時間(表2)，可以看出線性同時反演的計算效率遠遠高于鄰近算法-網格搜索方法。值得注意的是，隨著反演的微震事件個數的增加，線性同時反演方法所需的計算量幾乎成倍增長，而鄰近算法-網格搜索方法增長速度較慢。這是由于隨著微震個數的增加，線性同時反演需要求解的反演方程組也相應變大，而鄰近算法-網格搜索方法只在震源定位時增加了定位個數。因此，當實際數據中存在成百上千個微地震事件時，鄰近算法-網格搜索方法的適用性較強。

表2兩種方法的計算效率對比

4 結論

本文研究了微地震單井監(jiān)測情況下，利用鄰近算法-網格搜索方法反演震源位置和一維速度模型的可行性。主要討論了理論數據在有噪音和無噪音時的反演結果，以及搜索模型參數的數量對鄰近算法-網格搜索方法精度的影響。結果表明該方法具有較強的收斂性，但隨著模型參數的增加，反演誤差迅速上漲且計算量增加。為了驗證鄰近算法-網格搜索方法在實際數據中的應用效果，將反演結果與線性同時反演結果對比，結果表明線性同時反演對初始值的依賴性較強，容易陷入局部最優(yōu)解，而鄰近算法-網格搜索方法受初始值的影響較小，得到的震源定位結果分布在壓裂段附近，較為合理。

致謝：感謝東方地球物理公司提供的水力壓裂監(jiān)測數據以及審稿人提供的寶貴評閱意見。