地面與井中聯(lián)合觀測(cè)的微地震定位和速度結(jié)構(gòu)同時(shí)反演

黃國(guó)嬌　巴　晶　錢　衛(wèi)

(河海大學(xué)地球科學(xué)與工程學(xué)院地球探測(cè)研究所，江蘇南京 211100)

1　引言

水力壓裂技術(shù)是非常規(guī)油氣資源(如致密油氣、頁(yè)巖油氣)開(kāi)發(fā)過(guò)程中重要的增產(chǎn)措施，該技術(shù)通過(guò)將流體高壓注入地層，使巖石破裂以提高地層滲透率[1]。微地震監(jiān)測(cè)技術(shù)通過(guò)觀測(cè)、定位及分析微地震事件監(jiān)測(cè)生產(chǎn)活動(dòng)及地下?tīng)顟B(tài)，可用于水力壓裂效果分析[2-6]，在儲(chǔ)層開(kāi)發(fā)監(jiān)測(cè)中起到重要作用。目前，微震監(jiān)測(cè)技術(shù)應(yīng)用的主要困難在于： ①水力壓裂誘導(dǎo)的微地震事件通常信噪比較低，且觀測(cè)數(shù)據(jù)量大，初至旅行時(shí)難拾取； ②即使能夠較準(zhǔn)確地拾取初至旅行時(shí)(旅行時(shí)信息包含了地震事件的發(fā)震時(shí)刻和旅行時(shí)信息)，微地震事件的定位仍然存在誤差，原因是缺乏精確的地下速度模型。微地震反演中采用的速度模型一般根據(jù)地震、聲波測(cè)井、VSP等資料建立[6]，并采用壓裂前的射孔資料校正速度模型[7,8]。但是，射孔反演速度模型需已知射孔位置和發(fā)震時(shí)刻，而射孔事件的準(zhǔn)確發(fā)震時(shí)刻往往是未知的。再者，射孔數(shù)據(jù)一般在壓裂前采集，而在水力壓裂過(guò)程中，高壓流體會(huì)產(chǎn)生裂縫，新增裂縫和孔隙壓力的改變都會(huì)引起地層速度的變化[9]。因此如果采用射孔資料校正后的速度模型進(jìn)行微地震定位，可能存在較大誤差[10]。

針對(duì)微地震初至拾取困難的問(wèn)題，Song等[11]建立了一個(gè)微地震強(qiáng)信號(hào)的波形數(shù)據(jù)庫(kù)，通過(guò)分析強(qiáng)弱微地震信號(hào)的波形相關(guān)性提取弱微震事件的初至； 宋維琪等[12]提出了獨(dú)立分量分析與壓縮感知聯(lián)合理論的微地震弱信號(hào)的提取方法； 王晨龍等[13]和李振春等[14]采用地震干涉逆時(shí)定位算法，不需要拾取初至信息； 王鵬等[15]提出基于熵值的微震事件識(shí)別和恢復(fù)方法，以解決微地震監(jiān)測(cè)數(shù)據(jù)信噪比低的問(wèn)題； 李穩(wěn)等[16]將稀疏碼收縮法和小波變換相結(jié)合處理井下微地震數(shù)據(jù)，以保證有效地震信號(hào)的識(shí)別和提?。?曹俊海等[17]采用多道匹配跟蹤算法提高微地震信號(hào)的信噪比，進(jìn)而提高微地震信號(hào)識(shí)別的準(zhǔn)確性。

為得到較準(zhǔn)確的速度模型，減小速度模型對(duì)定位結(jié)果的影響，Grechka等[18]同時(shí)反演了介質(zhì)的各向異性參數(shù)和微地震位置，但假設(shè)地下為均勻介質(zhì)； Li等[19]將雙差定位法推廣至同時(shí)確定水平層狀各向異性參數(shù)及微地震定位； 譚玉陽(yáng)等[20]提出基于初至旅行時(shí)差的微地震速度模型反演方法，可在射孔起始時(shí)間未知的情況下反演地層速度模型。

在井下微地震監(jiān)測(cè)中，微地震事件與檢波器之間的距離通常在幾十至數(shù)百米的范圍內(nèi)，因此將地下介質(zhì)假設(shè)為水平層狀模型[7，21]，提出了地面與井中觀測(cè)旅行時(shí)聯(lián)合的微地震定位和速度結(jié)構(gòu)同時(shí)反演方法： 首先推導(dǎo)了旅行時(shí)關(guān)于速度、界面深度以及震源參數(shù)的一階偏導(dǎo)數(shù)，給出了Jacobi矩陣計(jì)算方法； 然后用最短路徑算法[22]追蹤射線路徑和計(jì)算旅行時(shí)；再采用共軛梯度算法求解帶約束的阻尼最小二乘問(wèn)題，更新地下層狀介質(zhì)的結(jié)構(gòu)(層速度和層界面深度)以及微地震定位(包括震源位置和發(fā)震時(shí)刻)； 最后討論了方法的抗噪聲能力。

2　方法原理

2.1　反演問(wèn)題

一般而言，基于旅行時(shí)數(shù)據(jù)的非線性同時(shí)反演問(wèn)題，可歸結(jié)為下列帶約束的阻尼最小二乘最優(yōu)化問(wèn)題，其反演問(wèn)題的L2范數(shù)目標(biāo)函數(shù)為

(1)

式中： δt=(δt1，δt2，…，δtM)是微地震事件的實(shí)際觀測(cè)旅行時(shí)與理論數(shù)據(jù)的旅行時(shí)殘差，M則為拾取的初至數(shù)量；δm=(λ1δV，λ2δD，λ3δX)，其中δV和δD分別是地下介質(zhì)各層速度和地層界面深度的相對(duì)改變量，δX為微地震事件的震源坐標(biāo)和發(fā)震時(shí)刻的相對(duì)改變量，λ1、λ2、λ3為各類參數(shù)的更新權(quán)重系數(shù)；A是Jacobi矩陣；μ是阻尼因子；a和b分別是反演解中更新各類參數(shù)更新量的上、下邊界值。

式(1)解的一般形式[23]為

[μCm+ATCdA]Zmδm=ATCdδt

(2)

式中：Cm、Cd為模型和數(shù)據(jù)空間的協(xié)方差矩陣；Zm為分段約束算子，即

(3)

式(2)是非線性問(wèn)題局部線性化的基本反演公式，其解具有局域解的特征。為了得到具有實(shí)際物理意義的解，可采用Cm、Cd及Zm的先驗(yàn)信息進(jìn)行約束。因此，選擇不同的組合(Cm、Cd、Zm)可得到不同形式的反演解。

本文采用Tarantola等[24]提出的廣義帶約束的阻尼最小二乘反演解，即Cm、Cd分別取模型和數(shù)據(jù)空間的協(xié)方差矩陣的逆，則解為

(4)

上述矩陣方程可基于共軛梯度法采用迭代算法求解[25]，其關(guān)鍵問(wèn)題是如何求解具有偏導(dǎo)性質(zhì)的Jacobi矩陣。

2.2　Jacobi矩陣元素計(jì)算

當(dāng)?shù)叵陆橘|(zhì)結(jié)構(gòu)(層速度和層深度)和微地震事件同時(shí)反演時(shí)，Jacobi偏導(dǎo)矩陣包含了三部分：一是旅行時(shí)關(guān)于速度變化的導(dǎo)數(shù)；二是旅行時(shí)關(guān)于層深度變化的導(dǎo)數(shù)； 三是旅行時(shí)關(guān)于微地震震源參數(shù)的導(dǎo)數(shù)。

(5)

式中： Δti,j代表第i個(gè)震源第j條射線的旅行時(shí)擾動(dòng)；vk和dk分別是第k層內(nèi)的速度分布和第k層界面的深度； Δvk和Δdk則分別是相應(yīng)速度和層界面深度的擾動(dòng)量；N為射線穿過(guò)的層數(shù)；K是與射線相交的界面數(shù)； Δxi,k(k=1，2，3，4) 是第i個(gè)微地震事件分別在x、y、z坐標(biāo)及發(fā)震時(shí)刻的擾動(dòng)量。

式(5)中第一項(xiàng)，即旅行時(shí)關(guān)于速度的一階偏導(dǎo)數(shù)為

(6)

式中L為射線在第k層內(nèi)的長(zhǎng)度。

式(5)中的第二項(xiàng)，即旅行時(shí)關(guān)于層界面深度變化的偏導(dǎo)數(shù)參照Rawlison[26]的方法根據(jù)圖1分析并推導(dǎo)。當(dāng)層界面的深度變化Δd時(shí)(從圖1中實(shí)線到虛線位置)，引起的旅行時(shí)變化為

(7)

則旅行時(shí)關(guān)于層界面深度變化的偏導(dǎo)數(shù)為

(8)

式中：向量wj、wj+1和wz如圖1所示；vj和vj+1則分別是層界面兩側(cè)的速度；e和f分別為界面擾動(dòng)后界面上、下兩側(cè)射線長(zhǎng)度的變化量。

式(5)中的第三項(xiàng)，即旅行時(shí)關(guān)于微地震事件參數(shù)(位置和發(fā)震時(shí)刻)變化的偏導(dǎo)數(shù)，可由以下步驟計(jì)算。首先，旅行時(shí)(t=ta-t0)關(guān)于發(fā)震時(shí)刻的導(dǎo)數(shù)為

(9)

其次，旅行時(shí)關(guān)于微地震事件位置變化的偏導(dǎo)數(shù)的計(jì)算，可參見(jiàn)圖2所示。若微地震震源的深度擾動(dòng)值為Δzh，則旅行時(shí)關(guān)于震源位置xh、yh、zh的偏導(dǎo)數(shù)[26]為

(10)

式中：θh如圖2所示；vh為震源處的速度值；φh、ψh分別為射線水平面投影后與x和y軸的夾角。

圖2　旅行時(shí)關(guān)于震源位置變化偏導(dǎo)數(shù)計(jì)算示意圖

3　數(shù)值模擬實(shí)例

為驗(yàn)證該算法的有效性及實(shí)用性，本文選用一個(gè)三維模型同時(shí)進(jìn)行地下介質(zhì)的層析成像及微地震事件的定位試驗(yàn)。如圖3所示，模型尺寸為200m×200m×50m，網(wǎng)格尺寸為4m×4m×5m。模型是四層的水平層狀介質(zhì)，速度從地表向下依次為2.0、2.4、2.8和3.2km/s，界面的深度分別為16、30和40m。為模擬地面和井中微地震，采用實(shí)際微地震中常用的地表“米”字型和井中聯(lián)合觀測(cè)方式(如圖3中三角形所示)，地表采集空間間隔為20m，井中采集空間間隔為10m，共46個(gè)檢波器。設(shè)定了4個(gè)虛擬微地震事件，其空間坐標(biāo)分別為S1(20m，40m，42m)、S2(100m，100m，45m)、S3(150m，180m，48m)、S4(170m，30m，26m)，激發(fā)時(shí)刻分別為10、15、20和5ms。

圖3　三維模型及微地震觀測(cè)系統(tǒng)三角形代表觀測(cè)點(diǎn)，圓點(diǎn)代表微地震事件

在理想條件下(正確的速度模型、無(wú)噪聲)采用該反演算法進(jìn)行微地震事件S1、S2和S3定位。定位結(jié)果如圖4所示，即使初始震源參數(shù)與真實(shí)值相差較大，反演的震源位置(圖4a～圖4c)和發(fā)震時(shí)刻(圖4d)也都與理論參數(shù)完全一致。

然而實(shí)際數(shù)據(jù)處理中不可能滿足理想條件，即地下介質(zhì)的速度結(jié)構(gòu)(層速度和層界面位置)未知、數(shù)據(jù)含噪聲等。下面給出本文反演算法同時(shí)進(jìn)行速度模型反演和微地震定位的結(jié)果，并討論其對(duì)噪聲的敏感性。首先選用兩種不同的微地震事件組合測(cè)試該反演算法同時(shí)進(jìn)行速度模型反演和微地震定位。第一種采用三個(gè)微地震事件(S1、S2和S3)的初至P波旅行時(shí)，第二種采用四個(gè)微地震事件(S1、S2、S3和S4)的初至P波旅行時(shí)，同時(shí)反演的結(jié)果如圖5(速度結(jié)構(gòu)反演結(jié)果)和圖6(微地震定位結(jié)果)所示。

圖4　理想條件下的微地震定位結(jié)果

從圖5 可以看出，當(dāng)僅采用三個(gè)微地震事件的旅行時(shí)數(shù)據(jù)時(shí)，第三層的速度結(jié)構(gòu)(層界面和層速度)層析的結(jié)果(圖5中藍(lán)線)與真實(shí)模型(圖5中紅線)差別較大；而當(dāng)采用四個(gè)微地震事件的旅行時(shí)數(shù)據(jù)時(shí)，層析結(jié)果明顯改善，無(wú)論是層速度值還是層界面位置(圖5中綠線)都接近真實(shí)值。出現(xiàn)這種情況是因?yàn)樵黾拥奈⒌卣鹗录4位于第二層介質(zhì)中，與微地震事件S1、S2和S3位置有較大區(qū)別，射線覆蓋及射線交錯(cuò)概率有所增加，提高了層析成像的分辨率。所以當(dāng)速度結(jié)構(gòu)反演結(jié)果得到改善時(shí)，相應(yīng)的微地震定位結(jié)果精度也有所提高，特別是震源深度(圖6b和6c)的定位結(jié)果。因此，當(dāng)射線的角度覆蓋較好時(shí)，同時(shí)進(jìn)行速度層析成像和微地震定位是可行的，反演精度較高。

圖5　不同微地震事件組合同時(shí)反演的速度結(jié)構(gòu)

圖6　不同微地震事件組合同時(shí)反演的微地震定位結(jié)果

考慮拾取初至旅行時(shí)可能存在誤差，在理論的初至P波旅行時(shí)數(shù)據(jù)中加入幅度為0.4ms的隨機(jī)誤差，四個(gè)事件同時(shí)反演的結(jié)果如圖7和圖8所示。結(jié)果顯示：加入隨機(jī)噪聲后，對(duì)微地震定位精度影響不大，且微地震震源水平位置的定位結(jié)果(圖8a)明顯好于震源深度的定位(圖8b和圖8c)。表1給出了具體的定位誤差，可見(jiàn)未加隨機(jī)噪聲之前水平定位誤差在幾厘米的范圍內(nèi)，加入隨機(jī)噪聲后水平定位誤差在十厘米左右，而深度的定位誤差為1～2m，但仍在可接受的誤差范圍內(nèi)。速度結(jié)構(gòu)層析的結(jié)果變化較大(圖7)，層速度和層界面深度都存在一定程度的偏差。這是由于當(dāng)速度結(jié)構(gòu)層析和微地震定位同時(shí)進(jìn)行時(shí)，兩者的反演結(jié)果實(shí)際上是一定程度的折中，并且層速度和層界面的深度在反演時(shí)也存在平衡的問(wèn)題?？傮w而言，微地震定位對(duì)旅行時(shí)中的隨機(jī)噪聲不敏感。

圖7　噪聲敏感性試驗(yàn)速度層析結(jié)果

x誤差/my誤差/mz誤差/mt0誤差/ms無(wú)噪加噪無(wú)噪加噪無(wú)噪加噪無(wú)噪加噪S10．002-0．013-0．026-0．1301．3991．975-0．077-0．109S2-0．0010．1680．0010．0500．8240．341-0．053-0．002S30．0140．118-0．0640．0240．8541．014-0．055-0．040S4-0．092-0．0620．1430．145-0．1690．233-0．103-0．227

圖8　噪聲敏感性試驗(yàn)微地震定位結(jié)果

上述的地面和井中聯(lián)合反演結(jié)果均迭代了20次左右。

4　結(jié)束語(yǔ)

為降低速度模型估計(jì)不準(zhǔn)對(duì)微地震定位精度的影響，本文提出了對(duì)水平層狀介質(zhì)同時(shí)進(jìn)行速度層析成像和微地震定位的反演算法。首先采用最短路徑算法計(jì)算初至波的射線路徑和旅行時(shí)，導(dǎo)出旅行時(shí)關(guān)于層速度、層界面深度及震源參數(shù)的偏導(dǎo)數(shù)，然后結(jié)合共軛梯度算法求解帶約束的阻尼最小二乘問(wèn)題，發(fā)展了一種同時(shí)反演速度結(jié)構(gòu)(層速度和層界面)和震源參數(shù)(震源位置和發(fā)震時(shí)刻)的方法。地面和井中聯(lián)合觀測(cè)條件下的數(shù)值模擬結(jié)果表明：該反演算法能有效反演地下介質(zhì)的速度結(jié)構(gòu)并進(jìn)行微地震定位，而且微地震定位結(jié)果對(duì)隨機(jī)噪聲不敏感。

本文基于模型正演無(wú)噪和加入隨機(jī)噪聲的數(shù)據(jù)，驗(yàn)證了所提出方法的有效性。在實(shí)際生產(chǎn)中，還存在地層起伏不平、介質(zhì)非均質(zhì)性強(qiáng)、數(shù)據(jù)信噪比低等諸多因素，因此還有必要結(jié)合實(shí)際數(shù)據(jù)與現(xiàn)場(chǎng)地質(zhì)情況，應(yīng)用、測(cè)試方法的可靠性和穩(wěn)健性，提高實(shí)際定位精度。

[1]張光明.水平井水力壓裂數(shù)值模擬研究[學(xué)位論文].安徽合肥：中國(guó)科學(xué)技術(shù)大學(xué)，2010.

Zhang Guangming.A Numerical Simulation Study on Hydraulic Fracturing of Horizontal Well[D].University of Science and Technology of China,Hefei,Anhui,2010.

[2]Rutledge J T and Phillips W S.Hydraulic stimulation of natural fractures as revealed by induced micro-earthquakes,Carthage Cotton Valley gas field,east Texas.Geophysics,2003,68(2)：441-452.

[3]Maxwell S C,Rutledge J,Jones R et al.Petroleum reservoir characterization using downhole microseismic monitoring.Geophysics,2010,75(5)：A129-A137.

[4]容嬌君,李彥鵬,徐剛等.微地震裂縫檢測(cè)技術(shù)應(yīng)用實(shí)例.石油地球物理勘探,2015,50(5)：919-924.

Rong Jiaojun,Li Yanpeng,Xu Gang et al.Fracture detection with microseismic.OGP,2015,50(5)：919-924.

[5]Li J L,Kuleli S,Zhang H J et al.Focal mechanism determination of induced microearthquakes in an oil field using full waveforms from shallow and deep seismic networks.Geophysics,2011,76(6)：WC87-WC101.

[6]張山,劉清林,趙群等.微地震監(jiān)測(cè)技術(shù)在油田開(kāi)發(fā)中的應(yīng)用.石油物探,2002,41(2)：226-231.

Zhang Shan,Liu Qinglin,Zhao Qun et al.Application of microseismic monitoring technology in development of oil field.GPP,2002,41(2)：226-231.

[7]Warpinski N R,Sullivan R B,Uhi J et al.Improved microseismic fracture mapping using perforation timing measurements for velocity calibration.Proceedings of the SPE Annual Technical Conference and Exhibition,Denver,Colorado,2003,84488.

[8]Van Dok R,Fuller B,Engelbrechet L et al.Seismic anisotropy in microseismic event location analysis.The Leading Edge,2011,30(7)：766-770.

[9]Block L V,Cheng C H,Fehler M C et al.Seismic imaging using microearthquakes induced by hydraulic fracturing.Geophysics,1994,59(1)：102-112.

[10]張曉林,張峰,李向陽(yáng)等.水力壓裂對(duì)速度場(chǎng)及微地震定位的影響.地球物理學(xué)報(bào),2013,56(10)：3552-3560.

Zhang Xiaolin,Zhang Feng,Li Xiangyang et al.The influence of hydraulic fracturing on velocity and microseismic location.Chinese Journal of Geophysics,2013,56(10)：3552-3560.

[11]Song F X,Kuleli H S,Toksoz M N et al.An improved method for hydrofracture-induced microseismic event detection and phase picking.Geophysics,2010,75(6)：A47-A52.

[12]宋維琪,李艷清,劉磊.獨(dú)立分量分析與壓縮感知微地震弱信號(hào)提取方法.石油地球物理勘探,2017,52(5)： 984-989,1041.

Song Weiqi,Li Yanqing,Liu Lei.Microseismic weak signal extraction based on the independent component analysis and compressive sensing.OGP,2017,52(5)：984-989,1041.

[13]王晨龍,程玖兵,尹陳等.地面與井中觀測(cè)條件下的微地震干涉逆時(shí)定位算法.地球物理學(xué)報(bào),2013,56(9)：3184-3196.

Wang Chenlong,Cheng Jiubing,Yin Chen et al.Microseismic events location of surface and borehole observation with reverse-time focusing using interfero-metry technique.Chinese Journal of Geophysics,2013,56(9)：3184-3196.

[14]李振春,盛冠群,王維波等.井地聯(lián)合觀測(cè)多分量微地震逆時(shí)干涉定位.石油地球物理勘探,2014,49(4)：661-666,671.

Li Zhenchun,Sheng Guanqun,Wang Weibo et al.Time-reverse microseismic hypocenter location with interferometric imaging condition based on surface and downhole multi-components.OGP,2014,49(4)：661-666,671.

[15]王鵬,常旭,王一博.基于時(shí)頻稀疏性分析法的低信噪比微地震事件識(shí)別與恢復(fù).地球物理學(xué)報(bào),2014,57(8)：2656-2665.

Wang Peng,Chang Xu,Wang Yibo.Automatic event detection and event recovery in low SNR microseismic signals based on time-frequency sparseness.Chinese Journal of Geophysics,2014,57(8)：2656-2665.

[16]李穩(wěn),劉伊克,劉保金.基于稀疏分布特征下的井下微地震信號(hào)識(shí)別與提取方法.地球物理學(xué)報(bào),2016,59(10)：3869-3882.

Li Wen,Liu Yike,Liu Baojin.Downhole microseismic signal recognition and extraction based on sparse distribution features.Chinese Journal of Geophysics,2016,59(10)：3869-3882.

[17]曹俊海,顧漢明,尚新民.基于局部相關(guān)譜約束的多道匹配追蹤算法識(shí)別微地震信號(hào).石油地球物理勘探,2017,52(4)：704-714.

Cao Junhai,Gu Hanming,Shang Xinmin.Microseismic signal identification with multichannel matching pursuit based on local coherence spectrum constraint.OGP,2017,52(4)：704-714.

[18]Grechka V,Singh P and Das I.Estimation of effective anisotropy simultaneously with locations of microseismic events.Geophysics,2011,76(6)：WC141-WC155.

[19]Li J L,Zhang H J,Rodi W L et al.Joint microseismic location and anisotropic tomography using differential arrival times and differential back azimuths.Geophy-sical Journal International,2013,195(3)：1917-1931.

[20]譚玉陽(yáng),何川,張洪亮.基于初至旅行時(shí)差的微地震速度模型反演.石油地球物理勘探,2015,50(1)：54-60.

Tan Yuyang,He Chuan,Zhang Hongliang.Microseismic velocity model inversion based on moveouts of first arrivals.OGP,2015,50(1)：54-60.

[21]Rutledge J T,Phillips W S and Schuessler B K.Reservoir characterization using oil-production-induced microseismicity,Clinton County,Kentucky.Tectonophysics,1998,289(1)：129-152.

[22]Bai C Y,Tang X P and Zhao R.2D/3D multiply transmitted,converted and reflected arrivals in complex layered media with the modified shortest path me-thod.Geophysical Journal International,2009,179(1)：201-214.

[23]Menke W.Geophysical Data Analysis：Discrete In-verse Theory.Academic Press Inc.,San Diego,California,1984,40-126.

[24]Tarantola A,Valette B.Generalized non-linear inverse problems solved using the least squares criterion.Reviews of Geophysics Space Physics,1982,20(2)：219-232.

[25]Zhou B,Sinadinovski C,Greenhalgh S A.Iterative algorithm for the damped minimum norm,least-squares and constrained problem in seismic tomography.Exploration Geophysics,1992,23(3)：497-505.

[26]Rawlison N and Sambridge M S.Seismic traveltime tomography of the crust and lithosphere.Advances in Geophysics,2003,46：81-197.