基于波動方程的微地震震源位置、震源時間及VTI介質各向異性參數聯合反演

薛清峰

1)中國科學院地質與地球物理研究所，油氣資源研究院重點實驗室，北京 100029 2)中國科學院地球科學研究院，北京 100029

0 引言

微地震監測技術作為一種有效監測水力壓裂情況以及儲層改造體情況的技術，近年來得到了飛速的發展(Thornton et al，2011)。微地震監測技術通過記錄致密儲層在水力壓裂過程中產生的地震響應來探知地下變形情況。對于微地震技術，其結果的準確度依賴于微地震震源位置與發震時間的準確度，而精確的微地震事件依賴于準確的背景物理場參數(Grechka et al，2017)。因此，如何獲得準確的微地震震源位置信息以及介質參數信息是水壓裂監測技術面臨的重要問題。

對于大多數的地震定位方法，識別地震震相以及拾取地震到時十分重要(Thurber et al，2000)。常規的微地震定位方法主要借助微地震信號的旅行時信息進行定位，例如，利用觀測數據與計算數據的旅行時差的方法(Luo et al，1991b)以及利用P波和S波到時差的方法(Lay et al，1995)。傳統的基于旅行時的方法需要對觀測數據進行有效地震事件的初至拾取，且對低信噪比的數據處理能力較差(Artman et al，2010)。近年來，人們借鑒地震勘探中的偏移成像原理，發展出不需要拾取震相走時信息且適用于低信噪比數據的類偏移定位方法(Artman et al，2010；Haldorsen et al，2013；O’Brien et al，2011；Sava，2011)。這類方法將微地震的震源類似成偏移成像中的繞射點，利用反射地震學中處理繞射點的成像方法進行震源的定位。相應的定位過程可以分為兩步，即首先延拓觀測數據“重構”地下波場，隨后施加合適的“成像條件”，得到震源位置和激發時間。逆時成像定位技術具有適應低信噪比數據、不需要拾取震相走時信息、較高的定位精度和可靠性等優點(Xue et al，2015、2016)。

對于非常規頁巖氣儲層，儲層各向異性發育明顯。因此，針對非常規油氣開發，不僅需要P波、S波速度場信息，還需要介質各向異性參數信息。當前速度和各向異性參數反演的前沿方法是全波形反演(FWI)方法，該方法通過目標函數的最小化來反演儲層介質參數模型。然而這樣的目標函數是高度非線性的，迭代往往會陷入局部極值中(Schuster，2017)。為了緩解這一問題，采取的方法是對數據進行框架式的表述，如采用初至旅行時時間(Luo et al，1991a)，可以反演獲得背景速度模型的中低波數至中間波數的細節(Sheng et al，2006)。骨架化的反演是準線性的，因此較傳統的全波形其具有更好的收斂性(Schuster，2017)。利用波動方程對旅行時間進行反演，被稱為波動方程旅行時反演(WT)方法(Luo et al，1991b)。微地震數據有2個主要特征，一是實際資料信噪比較低，二是獲得的數據是被動源數據，相較勘探地震的主動源數據，信號更加復雜(Berkhout et al，2011)。因此，利用其主要的骨架信息(旅行時間信息)進行背景介質參數場的反演，是更為合理可靠的方法。

震源參數的計算涉及介質參數信息，因此只有在獲得準確介質參數信息的情況下，才能得到精確的微地震震源參數。同時，準確的微地震參數信息也有助于介質參數場的反演。在地震學研究中，研究人員已經提出了基于射線方法的震源位置與介質參數聯合反演方法(Yuan et al，2016)，而利用更高精度的波動方程方法進行聯合反演還在起步階段(Zheng et al，2016)。因此，本文提出了基于波動方程的震源位置、震源時間及VTI介質各向異性參數聯合反演方法，利用該方法，VTI介質各向異性參數信息以及微地震事件的定位可以同時得到準確的求解。在微地震定位過程中，利用微地震介質參數反演得到介質參數模型，而在微地震介質參數反演過程中，所采用的震源參數信息則是由震源反演得到的，兩個過程相互促進，最終實現準確的求解。

本文首先介紹VTI介質的擬聲波方程，該方程為近似得到的解耦的偽聲波方程。在此基礎上，進一步提出基于VTI介質聲波方程的聯合反演理論。最后，利用本文所提的方法對合成資料進行試算，以驗證方法的有效性。

1 理論

VTI介質精確的頻散關系表達式為(Tsvankin，2001)

(1)

(2)

根據泰勒展開式的收斂條件，當P波和SV波頻散關系滿足式(3)時，式(2)所示的VTI介質下的純P波、純SV波方程有較好的近似

(3)

整理式(2)，解耦的VTI介質純P波方程在時間和波數域可表示為

(4)

其中，kx和kz為空間波數；p為時間和波數域的波場值。

在微地震事件的處理過程中，微地震事件的位置、事件的激發時間以及介質的模型參數均未知，并且這些參數之間存在高度的非線性關系。因此，如何有效地反演一個包含震源位置、震源激發時間、介質參數的聯合目標函數，并建立一套穩定且行之有效的反演方法至關重要。針對這一問題，本文給出的方法是建立一套如圖1 所示的聯合反演流程，即在反演框架內串行執行，首先進行震源位置的反演，隨后進行震源激發時間的反演，最后進行介質參數的反演，并將上述過程視為一次聯合反演階段，通過不斷迭代聯合反演階段，直至最終模型達到收斂要求或預期目標。

圖 1 聯合反演流程

1.1 基于逆時成像的微地震震源位置反演方法

借助波動方程逆時不變性的特點，將記錄波場作為邊界條件逆時延拓地震波場，使地震波場最終聚焦于震源處，實現對震源過程的成像(Steiner et al，2008)，這種確定震源參數的方式稱為震源逆時成像(Fink et al，2000)。以地震觀測數據作為輸入數據，通過波動方程逆時延拓地震波場，逆向重構地震波場在地下的傳播過程。二維情況下震源逆時成像的原理示意圖見圖2，圖中x-t平面內的雙曲線表示地震觀測記錄，紅色區域表示地震事件的空間位置與發震時刻，彩色曲面為地震波前擴展，顏色從藍色到紅色反映了能量的逐漸增強。

圖 2 震源逆時成像原理示意圖

在逆時重構過程中，波場會在發震時刻聚焦于震源位置，發震時刻過后又逐漸發散。因而，需要在波場延拓過程中施加合理的“成像條件”實現對震源成像，并獲取震源的相關屬性(如發震時刻、震源位置等)。通?？梢栽谝欢〞r間窗口范圍內選擇重構波場能量最大的點，將其視為推測的震源點(Gajewski et al，2005)。基于逆時成像的微地震震源位置反演方法可給出如下的震源位置目標函數

(5)

其中，p(t;xs，t0)表示(震源位置xs、震源激發時刻t0)產生的波場傳播到t時刻的波場值。

1.2 基于互相關的震源激發時間反演方法

激發時間更新主要考慮的是，獲得的微地震記錄長度可能與真正從微地震事件產生到檢波器接收的長度不一致，這就會對波動方程類的被動源成像方法造成誤差。其原因是在實際監測中，最終獲得的監測文件通常是單獨含有微地震事件的地震數據，而并非從開始施工監測直至接收到數據，因此，地震數據文件的開始記錄時間通常不能滿足逆時反傳成像的條件，即計算數據有可能無法反傳回發震時間零點，故需要借助互相關進行時間校正，在每次迭代中更新激發時間。更新激發時間的方法為，根據獲得的逆時成像定位結果I(x)，選取其中能量大于給定閾值的成像點，并將這些成像點推測為震源點，同時依據其能量的大小計算權重值；利用這些震源點以及權重值計算從震源點到接收點的波場，當檢波器得到重新計算的數據dcal后，利用式(6)的互相關方法，計算觀測數據dobs與計算數據dcal的互相關函數

(6)

其中，u(xr，t+τ，xs)obs、u(xr，t，xs)cal分別為觀測和計算的波場值。

使互相關函數f(xr，τ，xs)取得最大值的Δτ即為激發時間的更新量，從而得到最終更新后的激發時間Tcal=Tcal+Δτ。

由此，基于互相關的微地震震源激發時間反演方法可以給出震源激發時間目標函數

(7)

其中，p(xr，t+Δτ，xs)obs、p(xr，t，xs)cal分別為觀測和計算的波場值。

1.3 基于旅行時反演的介質各向異性參數反演方法

對于常規旅行時反演方法，通過引入各向異性參數VP0、ε和δ，可將其拓展到各向異性介質情況。介質參數反演則根據觀測數據，尋找使目標函數最小化的各向異性參數場。

旅行時反演中定義的目標函數為

(8)

其中，Δτ(xr，xs)=τobs(xr，xs)-τcal(xr，xs)，表示計算數據與觀測數據的初至旅行時差。若選擇最簡單的最速下降法來更新介質梯度，則第k次迭代過程中的介質參數更新公式可以表示為

m(x)(k+1)=m(x)(k)+αkγk(x)

(9)

其中，γk(x)表示參數更新的最速下降法方向，即負梯度方向；αk為步長，一般可通過線性搜索獲得。

目標函數S關于各向異性參數場m(x)的Fréchet導數為(Luo et al，1991b)

(10)

(11)

(12)

則式(11)中的分母可表示為

(13)

分子可以表示為

(14)

采用式(4)的VTI介質純P波方程進行波動方程旅行時反演，其3個各向異性參數VP0、ε和δ的梯度公式分別為(Feng et al，2016)

(15)

式中

(16)

由此，最終形成一個包含震源位置、震源激發時間、介質各向異性參數的聯合目標函數

(17)

2 數值算例

采用一個簡單的五層水平層狀模型(圖3)來驗證上述方法的有效性，該模型的水平方向有1000個網格點，垂直方向有1000個網格點，網格點間距為10m，共500個檢波器，均為地面接收，檢波器間隔20m，震源子波為50Hz雷克子波。真實的VTI各向異性參數模型如圖3(a)～(c)所示。為了更直觀地比較速度變換，提取X=600m位置處的一維各向異性真實模型與初始模型，如圖4 所示。

圖 3 五層水平層狀模型的各向異性參數(a)P波速度模型，紅色圓圈表示震源位置，粉色三角表示檢波器分布； (b)各向異性參數δ模型； (c)各向異性參數ε模型

圖 4 五層水平層狀模型各向異性真實模型與初始模型的比較(a)P波速度模型；(b)各向異性參數δ模型；(c)各向異性參數ε模型；黑色線代表各向異性真實模型；紅色線代表初始模型

由震源位置反演結果(圖5)可以看出，經過5次迭代后，定位結果與真實位置基本一致，證明了本文提出的聯合迭代反演定位方法的有效性。

圖 5 五層水平層狀模型震源位置反演結果(a)初始定位結果，紅色圓圈代表震源的真實位置，黃色十字代表能量聚焦的定位結果；(b)5次迭代的定位結果；色標表示無量綱相對能量值，顏色代表聚焦能量，值越大能量越強

震源時間反演結果隨迭代次數的變化情況，如圖6 所示，可以看到8個不同震源(編號S1～S8)的激發時間信息均得到了良好的反演，并且收斂速度較快。

圖 6 五層水平層狀模型震源激發時間反演結果藍色實線為反演結果；橙色虛線為真實值

3個各向異性參數的反演結果如圖7 所示，對比各向異性參數的真實模型(圖3)，可以看到本文所提的方法較好地還原了模型參數。同樣，為了更直觀地考察反演結果，提取了X=600m位置處的反演情況，如圖8 所示。從該一維視角可以看到，本文的反演結果較初始模型有了較大提升，細節刻畫更為清楚，每一層均得到較好的還原。

圖 7 五層水平層狀模型各向異性參數模型反演結果(a)參數VP0的結果；(b)參數δ的結果；(c)參數ε的結果

圖 8 五層水平層狀模型的各向異性真實模型、初始模型以及反演結果比較(a)P波速度模型；(b)各向異性參數δ模型；(c)各向異性參數ε模型；黑色線代表各向異性真實模型；紅色線代表初始模型；藍色線代表反演結果

由五層水平層狀模型介質參數反演的歸一化目標函數殘差(圖9)可以看到，收斂曲線在不斷下降，表明了聯合反演方法對VTI介質參數反演的有效性。

圖 9 五層水平層狀模型介質參數反演的歸一化目標函數殘差收斂曲線

進行簡單模型測試后，再將模型復雜化，利用一個推覆體模型進行計算。推覆體模型各向異性參數的真實模型如圖10(a)～10(c)所示，模型大小為：X方向800個網格點，Z方向180個網格點，空間網格10m；共200個檢波器，均為地表接收，40m間隔；8個震源，子波為50Hz雷克子波。

圖 10 復雜模型的各向異性參數(a)P波速度模型，黑色圓圈表示震源位置，粉色三角表示檢波器分布，震源S1～S8括號內的時間信息為震源激發時間； (b)各向異性參數δ模型；(c)各向異性參數ε模型

進行聯合反演的初始模型信息如圖11 所示。經過20次迭代反演，反演定位結果(圖12(b))與真實位置基本一致，進一步驗證了本文的聯合反演方法在提高微地震事件定位結果精度上的有效性與可靠性。

圖 11 復雜模型的各向異性參數反演初始模型(a)P波速度初始模型；(b)各向異性參數δ初始模型；(c)各向異性參數ε初始模型

圖 12 復雜模型的震源位置反演結果(a)初始定位結果，紅色圓圈代表震源的真實位置，黃色十字代表能量聚焦的定位結果；(b)20次迭代的定位結果；色標表示無量綱相對能量值，顏色代表聚焦能量，值越大能量越強

8個不同震源(編號S1～S8)的時間反演結果如圖13 所示，可以看到，在復雜模型且多震源不同激發時間的情況下，本文的聯合反演方法同樣較快速地收斂到了真實的激發時間。

圖 13 復雜模型的震源激發時間反演結果藍色線為反演結果；橙色線為真實激發時間

復雜模型下3個各向異性參數的反演結果如圖14 所示，對比真實的各向異性參數(圖10)可見，本文的方法對VTI各向異性參數信息有較好的反演。由復雜模型介質參數反演的歸一化目標函數殘差(圖15)可見，收斂曲線不斷下降，進一步說明了本文的聯合反演方法在介質參數反演中的有效性。

圖 14 復雜模型的各向異性參數模型反演結果(a)參數VP0 的結果；(b)參數δ的結果；(c)參數ε的結果

圖 15 復雜模型介質參數反演的歸一化目標函數殘差收斂曲線

3 結論

對于非常規油氣開發，準確的微地震震源位置信息是關乎水力壓裂施工的重要支撐數據，而微地震震源位置的精確提取依賴于準確的激發時間和可靠的介質參數。在波動方程逆時成像定位以及波動方程旅行時反演的基礎上，本文提出了同時反演震源信息(位置、時間)與介質參數信息的聯合反演方法，可以獲得精確的微地震震源信息以及介質信息，并通過數值模型檢驗了該方法的有效性和穩定性。本文建立的微地震聯合反演方法計算效率較高，易于實現，在實際應用中可對微地震的震源信息的提取及介質各向異性模型的建立提供有效信息。