基于波動(dòng)方程的微地震震源位置、震源時(shí)間及VTI介質(zhì)各向異性參數(shù)聯(lián)合反演

薛清峰

1)中國(guó)科學(xué)院地質(zhì)與地球物理研究所，油氣資源研究院重點(diǎn)實(shí)驗(yàn)室，北京 100029 2)中國(guó)科學(xué)院地球科學(xué)研究院，北京 100029

0 引言

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

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

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

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

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

1 理論

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

(1)

(2)

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

(3)

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

(4)

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

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

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

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

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

圖 2 震源逆時(shí)成像原理示意圖

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

(5)

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

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

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

(6)

其中，u(xr，t+τ，xs)obs、u(xr，t，xs)cal分別為觀測(cè)和計(jì)算的波場(chǎng)值。

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

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

(7)

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

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

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

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

(8)

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

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

(9)

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

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

(10)

(11)

(12)

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

(13)

分子可以表示為

(14)

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

(15)

式中

(16)

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

(17)

2 數(shù)值算例

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖 12 復(fù)雜模型的震源位置反演結(jié)果(a)初始定位結(jié)果，紅色圓圈代表震源的真實(shí)位置，黃色十字代表能量聚焦的定位結(jié)果；(b)20次迭代的定位結(jié)果；色標(biāo)表示無(wú)量綱相對(duì)能量值，顏色代表聚焦能量，值越大能量越強(qiáng)

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

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

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

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

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

3 結(jié)論

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