基于微分格式的微地震走時反演方法研究

胡 飛， 蘇 奧

(1.中國科學技術大學 地球與空間科學學院，合肥 230026；2.東方地球物理公司 研究院，涿州 072750)

?

基于微分格式的微地震走時反演方法研究

胡 飛1，2， 蘇 奧2

(1.中國科學技術大學 地球與空間科學學院，合肥 230026；2.東方地球物理公司 研究院，涿州 072750)

在頁巖氣開發涉及的微地震監測中，只利用聲波測井資料建立的速度模型并不足夠準確，通常的做法是利用測井資料與射孔數據、井下爆炸索數據、井中下落球震數據結合的方法構建模型，但射孔、爆炸索等數據的震源初始時間很難準確測量。初始時間的準確性將會影響速度模型的準確性，進一步影響定位結果的準確性。這里提出一種不需要震源初始時間的反演方法，與傳統的速度反演方法相比，該方法基于不同震相的到時差信息。對同一震源，提取P-P，S-S和P-S震相到時差，這三部分基于微分格式的到時差信息，都可完全消除初始時間的影響。對合成數據測試、分析和討論的結果，證實了該方法實現速度模型校正的可行性，并給出了該方法在實際數據應用中的示例。

微地震； 微分格式； 走時反演； 速度模型校正

0 引言

微地震監測技術可以追溯到上世紀70年代，在過去幾十年里已經被證明是一種了解地下結構等信息的不可或缺的技術。在微地震成像中，速度模型的獲取最為重要，因為更加準確的目標區域速度結構，可以提高微地震事件定位的精確度。另外在微地震定位結果的不確定性分析中，速度模型誤差帶來的定位錯誤要比走時拾取誤差帶來的定位錯誤大得多[1]。考慮到頁巖氣所在區域的沉積地質特性，工業中通用的做法是利用聲波測井資料建立層狀的速度模型。但是這也帶來一個重要問題——該速度模型并不是足夠準確的，這將會影響到后續處理，比如地震定位。速度模型不夠準確主要有四方面原因，①頻率不同，聲波測井的頻率可以高達幾千赫茲，而微地震的頻率一般只有幾十赫茲，更高的頻率帶來更大的衰減系數，導致測量不準；②利用測井資料獲取的速度信息只是縱向的速度，實際需要的是穿過目標區域的橫向速度信息，目標區域的橫向速度值相對于測井資料結果，有10%～20%的速度誤差[2]；③射孔數據，井下爆炸索數據的初始時間很難準確測量[2-3]，Warpinski等[4]研究發現是由于安全射孔工藝技術導致的，該技術中的電容器放電和火花間隙會導致導火索點火被延遲；④水力壓裂過程會改變目標區域的速度結構，比如該區域已經進行過一段時間的生產，會導致壓力下降，或者地下流體改變，這些都會改變該區域的速度，測井得到的資料相對于生產監測時，速度信息已經不再正確[5]。因此在微地震監測工作開始前，必須要進行速度模型的校正。

一些評估水力壓裂速度模型的方法已經被提出，Block等[6]提出一種參數分離方法來聯合反演破裂區速度信息、震源位置、發震時刻；Warpinsk等[4]基于射孔數據和初始時間的評估，提出一種最小二乘的線性反演方法。這些方法都提高了微地震成像的效果，但是他們的方法依然存在初始時間無法精確測量的問題。Pei等[5]利用快速模擬退火法和線性奧卡姆反演方法去校正速度模型。他們把P、S震相利用聯合目標函數結合起來，不過初始時間無法準確測量的問題依舊存在。

這里提出一種基于微分格式的非線性最小二乘反演方法，通過擬合微分格式的旅行時而不是直接擬合旅行時，來校正測井資料構建的速度模型。利用微分格式的走時信息，其優勢在于能很好解決震源初始時刻難以準確測量的問題。將該方法應用到合成數據和實際數據的實例中，可以實現：只反演P波速度、只反演S波速度、同時反演P波和S波速度。

1 反演原理和方法

在頁巖氣的工業生產中，由于研究尺度較小，地層多為沉積特性，故通常建立層狀結構模型。該模型多為結合測井資料與射孔，爆炸索等數據獲得，但是存在震源起爆初始時間很難精確測量的問題。基于微分格式的微地震走時反演方法，試圖解決這一難題。該方法通過建立最小二乘目標函數，既實現對微分格式走時信息的擬合，也實現對模型平滑的約束。以測井資料提供的信息為初始速度結構，迭代更新速度模型，至目標函數前、后兩次殘差達到最小且保持不變時停止迭代，從而達到校正測井速度模型的目的。

1.1 目標函數

通過擬合檢波器記錄的地震波走時和速度模型正演計算的走時，我們可以獲取地下結構的層速度，作者提出的反演方法是通過擬合微分格式的走時而不是直接擬合走時來校正速度模型。具體來講，對同一震源，考慮不同檢波器P-P震相到時差，S-S震相到時差，P-S震相到時差，這三部分基于微分格式的到時差信息，都可完全消除初始時間的影響，所以目標函數中要加入這三個震相時差的約束項。考慮到測井資料模型雖然不夠精確，但已是很好的初始速度模型，也是較為準確的真實速度結構，所以目標函數要加入初始速度模型(測井資料)的約束項。同時考慮到問題的穩定性和模型的平滑性，我們應用吉洪諾夫正則化[7]到目標函數中，并選擇L2范數作為最小二乘評估算子。綜上所述，可采取如下具有物理意義的目標函數Φ(m)：

Φ(m)=a0‖(Δdp-ΔG(mp))‖2+a1‖(Δds-ΔG(ms))‖2+a2‖[(ds-dp)-(G(ms)-G(mp))]‖2+τp‖L(mp-mp0)‖2+τs‖L(ms-ms0)‖2

(1)

式中：Δdp和ΔG(mp)分別是觀測和計算的P-P震相走時微分；Δds和ΔG(ms)分別是觀測和計算的S-S震相走時微分；mp和ms分別是P波和S波慢度；dp和ds分別是觀測的P波到時、S波到時，都包括了震源初始時間；G(mp)和G(ms)分別是利用P波慢度和S波慢度計算的合成到時；a0,、a1和a2分別是P-P、S-S、P-S震相權重因子；τp和τs是模型平滑控制參數；mp0和ms0分別是從聲波測井數據獲取的初始P波慢度和初始S波慢度；L是正則化算子。

目標函數Φ中前三項分別給出了P-P、S-S和P-S震相的走時殘差信息(觀測到的地震波震相走時差與計算得到的震相走時差的殘差)，后兩項分別給出了對P、S波的新速度模型相對測井提供的初始模型變化量的約束。算子L采用二階Laplace算子，考慮到模型空間或者實際地球介質的速度結構存在非線性變化，所以采用具有非線性插值約束的二階梯度Laplace算子是合理的。通常來講，平滑控制因子τp和τs的選取具有一定的主觀性，初始值可以給1.0。權重因子a0、a1和a2可以根據地震資料各震相的信噪比選擇，或者根據用戶擬合某些震相的要求選取權重大小。在此作者根據P-P、S-S和P-S三震相的走時殘差均方根不同，提供一種求解方案，具體公式如下：

(2)

即a0=1，a1取P-P震相與S-S震相旅行時差的均方根之比，a2取P-P震相與P-S震相旅行時差的均方根之比，其中RMS代表求解均方根值。這樣選擇權重的出發點是三組微分格式走時中，P-S震相到時差的量級要比P-P震相到時差，S-S震相到時差大很多。

1.2 反演方法

求解使得目標函數取值最小的模型值即新的速度結構，我們使用最小二乘的分析方法來求解目標函數的極值問題，這需要進行兩部分的工作：①正演計算炮點到檢波點的走時；②求解目標函數的極值。

在第一部分正演問題中需要計算炮點到檢波器的走時，傳統上有兩種方案計算旅行時:射線追蹤法和波前追蹤法。射線追蹤包括試射法和彎曲法，波前追蹤包括有限差分法和圖法。試射法和彎曲法是追蹤射線從源點到檢波器點路徑的方法，有限差分是用程函方程去求解走時的方法，圖法是基于惠更斯原理的方法。在微地震合成數據的測試中，模擬微地震生產中產生射孔數據的觀測系統，首先建立層狀模型，在井下進行震源激發，在井中布置檢波器接收，并利用圖法中的最短射線路徑法(SPR)[8]計算P波、S波旅行時。考慮到實際資料處理是從檢波器記錄的波形上拾取到時信息作為觀測走時，所以為了結論的嚴謹性和正確性，在合成數據測試中，采用廣義反透射系數方法(GRTM)[9]先計算波形，然后手動拾取初至走時作為觀測數據，在迭代過程中使用SPR方法得到計算走時數據。

在第二部分求解目標函數極值的過程中，為避免計算繁雜的海賽矩陣，通過Gauss-Newton (GN) 法對目標函數進行線性化處理，以期求得迭代形式的更新公式。首先對目標函數(公式)進行形式變換，利用L2范數的定義，目標函數可改寫為：

(3)

其中：M是檢波器記錄的總道數；N是模型層數。公式中的各矩陣定義如下：

公式(3)相對公式(1)形式上得到大幅簡化，很好地將觀測數據d與計算數據G進行了結合，這將使得反演過程得到簡化。根據GN方法，當目標函數Φ(m)取得極值時，有遞推公式如下：

(4)

對公式(3)求導，可求得如下結果：

(5)

將公式(5)代入公式(4)可得到遞推公式(6)。

(6)

在公式(6)中，各矩陣定義如下：

在矩陣Ak中,第i行第j列元素代表了第i條射線在模型第j層中的射線長度。矩陣Wk是對不同數量級的P-P、S-S和P-S震相數據進行權重控制的因子。考慮到問題的非線性情況和求解的穩定性，我們在GN方法得到的公式(6)左邊加入可調節的阻尼項εkI。對于阻尼最小二乘法，最大的困難是如何選取阻尼εk值。麻省理工學院的Morgan教授[10]提出了一種方法:εk=α×σrms， 其中α是一個經驗性的參數(大約0.1)，σrms是公式(6)右邊部分的均方根誤差。這個方法的優勢在于隨著每次迭代的進行，阻尼因子εk會一直變化并且被自動賦值。具體來講，如果目標函數的值不是足夠小或者擬合效果較差時，σrms將會變大，這進一步導致了一個足夠大的εk被自動應用到阻尼系數中，最終會導致相對小的模型改變。另一方面，隨著反演過程進行，結果越來越收斂，σrms就會減小，即一個足夠小的εk將會被自動應用到阻尼系數中。該方法在一些已經發表的文章中已經被成功應用[11-13]。綜上所述我們得到最終的迭代公式：

(7)

我們可以應用共軛梯度法等方法通過多次迭代來更新速度模型，迭代過程中求解觀測數據和計算數據殘差的均方根誤差，當該均方根誤差降低到穩定不變時迭代結束，迭代所得結果即新速度模型。

基于微分格式的走時反演可以分為三種情況，①只反演P波速度；②只反演S波速度；③同時反演P波和S波速度。在上面的公式推導中我們以同時反演P波和S波為例給出了公式(1)，對于另外兩種情況，我們只需將對應的系數和權重因子設為“0”。

2 合成數據測試

通過上面的理論分析，我們將用合成數據進一步闡釋該方法的可行性。圖1是在文獻[14]中使用的德克薩斯州微地震工區模型基礎上，建立的一個微地震層狀模型。觀測系統由2個射孔震源和24個檢波器組成(圖1(a))，其中檢波器序列位于垂直井中，間隔15 m，深度從2 070 m至2 415 m。圖1(b)展示了震源到檢波器的射線路徑,圖中震源位置為(400 m,2 600 m)和(800 m,2 600 m)。

表1給出了該速度模型的具體層數，厚度，P波和S波真實速度, 其中3 000 m以下為無限空間。從表1中可以看到，該模型有高速層，有低速層，在目標區域(可能存在頁巖氣的區域)附近的層位中(深度2 171 m～2 457 m)層厚最小34 m，最大119 m，目標區域厚度大約300 m，符合常見頁巖氣地層的主要速度特征。

測井數據提供的模型很接近真實結構，但通常與目標區域橫向速度存在10%～20%[2]的誤差。因此在合成數據測試中，將真實模型每層的速度進行±[10%, 20%]隨機擾動后，可當作測井資料提供的速度信息，同時也可作為反演的初始模型。然后利用微分格式的走時反演方法，進行只反演P，只反演S，同時反演P、S速度。圖2顯示利用GRTM方法計算的兩個震源A和B的波形，紅色豎線和藍色豎線分別是手動拾取的P波、S波到時。圖2中橫軸代表走時，單位s,縱軸是檢波器序號，每個檢波器記錄中的三種顏色，代表了三分量的波形記錄。

圖1 速度速度信息及射線路徑

圖3分別給出了只反演P、只反演S和同時反演P、S速度的結果。圖3中橫軸是速度，縱軸是深度，虛線是初始速度模型，紅線是真實速度模型，藍線是反演結果。從圖3中可以看到，三種反演方案在目標區域層都得到了很好的速度恢復，而頂、底兩層的恢復效果較差，從圖1(b)看到射線路徑未經過頂、底兩層，故無法進行這兩層的速度校正，這與圖3中結果一致。而圖3顯示頂、底兩層的速度仍發生了變化，這是因為目標函數中二階Laplace算子對模型慢度參數的運算所致。在實際處理中，考慮到該運算對慢度改變量的未知性，建議這類層位在校正后仍沿用測井資料提供的速度信息。

表1 模型參數

圖2 利用GRTM方法計算得到兩個震源的波形

圖3 不同波形反演結果

圖4 震源初始時間T0的影響

3 分析與討論

基于微分格式的微地震走時反演方法主要特點，解決震源起爆初始時間難以準確測量的困難，以測井資料提供的初始速度結構為約束，來進行地震波層速度的校正。考慮到初始時間和初始模型在反演過程中的重要性，將對初始時間T0和初始速度模型對反演結果的影響分別討論。

3.1 震源初始時間T0對反演結果的影響

微地震研究區域通常較小，同時在井下激發和接收，使得微地震觀測走時量級上偏小，較小的初始時間誤差也會對速度模型的正確性產生較大影響。圖4給出了震源初始時間T0對反演結果的影響。測試中，對正確T0加上[-2,2] ms隨機擾動作為錯誤T0。圖4的結果顯示，較小的震源初始時間偏差，就會導致利用絕對走時反演得到的速度模型出現錯誤，這進一步說明了初始時間的重要性，同時可以看到作者提出的方法在震源初始時間存在錯誤時，仍然可以得到準確的速度結構，表明了該方法校正速度模型的可行性。

3.2 初始模型對反演結果的影響

在合成數據的測試中，用來模擬測井資料的初始速度結構是通過真實模型隨機擾動得到的。為了進一步說明反演方法的可行性，進行了不同常速度初始模型的討論。圖5給出了不同常速度初始模型所對應的反演結果。從圖5的結果可以看到，基于微分格式的反演方法對初始模型并沒有太強的依賴，對初始模型有一定的抵抗性。圖5中初始速度模型分別為500 m/s、 1 000 m/s、4 000 m/s、8 000 m/s。

4 應用實例

將微分格式走時反演方法應用到一個微地震的下落球震(dropball data)實際資料中，并給出相應的反演結果。圖6是從聲波測井資料得到的初始速度模型(速度-深度圖)，圖6中左邊為S波速度曲線，右邊為P波速度曲線。圖7(a)給出了此次工程中四個球震震源和12個檢波器組成的觀測系統圖，

圖5 初始速度模型對反演結果影響

圖7中星號(*)為四個下落小球震源，倒三角(▽)是12個檢波器。由于是層狀模型，將3維震源坐標 旋轉到x-z平面內后，四個震源坐標從左至右依次為(118.56 m，1 423.44 m)，(230.29 m,1 419.96m)，(534.41 m,1 422.93 m)，(863.86 m,1 425.73 m)。12個檢波器坐標從1 120 m至1 285 m，間隔為15 m。

圖6 由聲波測井得到的P，S波初始速度模型

此次研究的目標區域位于深度1 100 m至1 450 m，將圖7(a)中目標區域(紅色方框所示)放大得到圖7(b)，利用SPR計算的射線路徑也展示在圖7(b)中。目標區域以外并無射線經過，所以建議這些無射線經過的層位最終速度仍沿用測井資料提供的速度信息。

圖7 實際資料模型和射線路徑

圖8顯示了12個三分量檢波器記錄到的波形信息，這里只展示檢波器記錄到的最左邊震源的波形。圖8中橫軸代表走時，單位s，縱軸是檢波器序號，每個檢波器記錄中的三種顏色，代表了三分量的波形記錄，紅豎線和藍豎線分別為手動拾取的P、S波初至時間。

圖8 球震數據實際波形和手動拾取的P、S到時

圖9給出了實際數據只反演P、只反演S、同時反演P、S的結果。由于不知道真實速度模型，故圖9中只有初始模型和反演結果。

對于實際數據，由于不知道真實的地下結構，雖然不能通過反演后走時曲線的擬合圖來判斷速度結構正確性，但是走時擬合圖可以作為檢驗結果正確性的一種輔助方式。圖10 (分別給出了只反演P、只反演S、同時反演P和S時走時曲線擬合圖，在圖中上下左右四塊分別對應了沿x方向的四個震源。考慮到實際數據處理中存在噪音，拾取誤差等，可認為走時曲線擬合效果較好。

5 結論

這里介紹了頁巖氣開采中，以測井資料為基礎，校正微地震速度模型的一種可行方案。利用微分格式的走時信息，很好地解決了射孔數據等震源初始時間不能測量或者很難準確測量的困難。作者盡可能組合更多走時信息：P-P震相到時差、S-S震相到時差、P-S震相到時差，這三部分信息都可以排除初始時間的影響。合成數據的測試結果，以及初始時間、初始模型的討論結果，都證明了該方法的可行性。還通過一個球震實際資料給出了實際應用的示范。該方法主要有兩方面的貢獻：①解決了震源初始發震時間很難準確測量的問題，②將旅行時以微分格式應用到速度反演中。

圖9 實際球震數據反演結果圖

圖10 四個震源的走時曲線擬合圖

[1] JONES G A, KENDALL J M, BASTOW I D, et al. Locating microseismic events using borehole data [J]. Geophysical Prospecting, 2014, 62(1): 34-49.

[2] WARPINSKI N. Microseismic monitoring: Inside and out [J]. Journal of Petroleum Technology, 2009, 61(11): 80-85.

[3] PEI D, QUIREIN J A, COMISH B E, et al. Velocity calibration for microseismic monitoring: Applying smooth layered models with and without perforation timing measurements [C]. SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers, 2008.

[4] WARPINSKI N R, SULLIVAN R B, UHL J, et al. Improved microseismic fracture mapping using perforation timing measurements for velocity calibration [J]. SPE Journal, 2005, 10(01): 14-23.

[5] PEI D, QUIREIN J A, COMISH B E, et al. Velocity calibration for microseismic monitoring: A very fast simulated annealing (VFSA) approach for joint-objective optimization [J]. Geophysics, 2009, 74(6): WCB47-WCB55.

[6] BLOCK L V, CHENG C H, FEHLER M C, et al. Seismic imaging using microearthquakes induced by hydraulic fracturing [J]. Geophysics, 1994, 59(1): 102-112.

[7] TIKHONOV A N, ARSENIN V I A, JOHN F. Solutions of ill-posed problems [M]. Washington, DC: Winston, 1977.

[8] MOSER T J. Shortest path calculation of seismic rays [J]. Geophysics, 1991, 56(1): 59-67.

[9] CHEN X. Seismogram synthesis for multi- layered media with irregular interfaces by global generalized reflection/transmission matrices method. I. Theory of two-dimensional SH case [J]. Bulletin of the Seismological Society of America, 1990, 80(6A): 1696-1724.

[10] MORGAN F D O. Electronics of sulfide minerals: implications for induced polarization [D].Massachusetts Institute of Technology, 1981.

[11] MACKIE R L, MADDEN T R. Three- dimensional magnetotelluric inversion using conjugate gradients [J]. Geophysical Journal International, 1993, 115(1): 215-229.

[12] ZHANG J, MACKIE R L, MADDEN T R. 3-D resistivity forward modeling and inversion using conjugate gradients [J]. Geophysics, 1995, 60(5): 1313-1325.

[13] ZHANG J, TOKSOZ M N. Nonlinear refraction traveltime tomography [J]. Geophysics, 1998, 63(5): 1726-1737.

[14] WONG J, MANNING P M, HAN L, et al. Synthetic microseismic datasets . Focus, 2011.

Microseismic travel-time inversion based on time difference

HU Fei1,2, SU Ao2

(1.School of Earth and Space Science, University of Science and Technology of China, Hefei 230026, China;Geophysical Research Institute, BGP, CNPC, Zhuozhou 072750, China)

During microseismic monitoring for unconventional shale gas exploitation, initial velocity model is constantly built from sonic data or a combination of sonic data and perforation shot (or string shot, drop ball). Unfortunately, it's quite difficult to obtain the origin time. The origin time of perforation shot affects the accuracy of velocity structure, consequently the accuracy of subsequent location results. Therefore, we demonstrate an inversion procedure to calibrate the velocity model without origin time but using differential arrival times. Comparing with traditional methods, we combine travel-time data, such as differential P-P arrival times, differential S-S arrival times, and differential S-P arrival times. Both synthetic examples and real data results show the improvement and effectiveness of the inversion method in velocity model calibration.

microseimic; differential; travel-time inversion; velocity calibration

2014-07-25 改回日期：2014-10-11

胡飛(1988-)，男，碩士，從事微地震地球物理研究，E-mail:hufei@mail.ustc.edu.cn 。

1001-1749(2015)04-0478-10

P 631.4

A

10.3969/j.issn.1001-1749.2015.04.11