薛清峰

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

0 引言

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

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

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

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

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

1 理論

VTI介質(zhì)精確的頻散關(guān)系表達式為(Tsvankin，2001)

(1)

(2)

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

(3)

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

(4)

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

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

圖 1 聯(lián)合反演流程

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

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

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

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

(5)

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

1.2 基于互相關(guān)的震源激發(fā)時間反演方法

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

(6)

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

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

由此，基于互相關(guān)的微地震震源激發(fā)時間反演方法可以給出震源激發(fā)時間目標函數(shù)

(7)

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

1.3 基于旅行時反演的介質(zhì)各向異性參數(shù)反演方法

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

旅行時反演中定義的目標函數(shù)為

(8)

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

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

(9)

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

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

(10)

(11)

(12)

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

(13)

分子可以表示為

(14)

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

(15)

式中

(16)

由此，最終形成一個包含震源位置、震源激發(fā)時間、介質(zhì)各向異性參數(shù)的聯(lián)合目標函數(shù)

(17)

2 數(shù)值算例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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