宋維琪,徐奔奔,楊勤勇,郭全仕,姜宇東,喻志超,秦 晅

(1.中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島266580;2.中國石油化工股份有限公司石油物探技術(shù)研究院,江蘇南京211103)

震源反演問題源于天然地震。微地震震源反演問題和天然地震震源反演問題本質(zhì)上是相似的,但由于形成震源的機(jī)制不同,因此其具體實(shí)現(xiàn)方法不同。天然地震震源是構(gòu)造活動引起地層錯(cuò)動而產(chǎn)生強(qiáng)烈的地震波,屬于斷裂活動震源。微地震震源是壓裂后巖石破裂,釋放能量產(chǎn)生地震波[1]。在接收點(diǎn)觀測的微地震記錄是地層、震源和其它干擾因素的綜合響應(yīng)疊加結(jié)果。一般在地層速度均勻或橫向均勻、震源為均勻輻射點(diǎn)源(輻射能量大小在各個(gè)方向相同)情況下,利用觀測記錄計(jì)算的事件的方位角和震源的方位角相同。如果地層速度橫向變化或震源為非均勻輻射,則事件的方位角和震源的方位角不同。要實(shí)現(xiàn)事件的高精度定位,必須去除或校正掉地層的橫向變化效應(yīng)和震源輻射的非均勻效應(yīng)。研究源輻射的非均勻性問題,屬于震源反演問題[2]。

正演是反演的基礎(chǔ),考慮到實(shí)際反演需要,一般采用射線追蹤方法進(jìn)行正演。對于寬入射角接收的微地震射線追蹤正演,基本采用走時(shí)模擬,基于射線追蹤[3]的微地震高精度振幅、相位正演模擬很少,而這直接又關(guān)系到方位角反演的準(zhǔn)確性?？碧降卣痤I(lǐng)域研究的AVO反射系數(shù)近似計(jì)算公式,是在限定角度下推導(dǎo)的。對于寬入射角透射、反射系數(shù)的計(jì)算,還沒有實(shí)際有效的方法。本文采用簡化的高斯束射線追蹤方法[4]進(jìn)行正演。

隨著微地震監(jiān)測技術(shù)的發(fā)展,關(guān)于微地震震源機(jī)制及反演問題的研究也逐漸增多,但大多只注重某一方面。本文從地層響應(yīng)效應(yīng)、檢波器定位及震源等方面進(jìn)行系統(tǒng)研究,提出了新的研究思路和技術(shù)方法。

1 地震震源矩張量、張量分析

矩張量反演已成為現(xiàn)今天然地震主要的震源機(jī)制求解方法[5],一般采用力偶源組合震源模型。矩張量存在9個(gè)矩張量分量,表示為

(1)

矩張量具有對稱性,因此存在6個(gè)獨(dú)立分量。如果能夠獲得這6個(gè)獨(dú)立分量的值,就可以確定震源機(jī)制的3個(gè)軸(P,T,B軸)；同時(shí),根據(jù)震源機(jī)制軸面之間的位置關(guān)系,相應(yīng)的節(jié)面空間位置也得到確定。矩張量反演就是為了確定6個(gè)獨(dú)立矩張量的值,震源傳播方程可以寫作[6]

(2)

其中,U是地震波傳播到地面時(shí)檢波器接收到的位移序列;格林函數(shù)G描述地層介質(zhì)因素;M是6個(gè)獨(dú)立分量序列。(2)式把記錄到的地震信號很好地分成了震源和地層介質(zhì)兩個(gè)因素。通過記錄到的位移和格林函數(shù)可以反演得到6個(gè)獨(dú)立的矩張量。6個(gè)獨(dú)立的矩張量如果為純雙力偶的情況,就可以求解矩張量的3個(gè)特征向量,這3個(gè)特征向量對應(yīng)震源機(jī)制的P,T,B3個(gè)軸。

震源可以看成是一個(gè)不連續(xù)的點(diǎn),也可以看作一個(gè)雙力偶源。作為一個(gè)雙力偶源時(shí),可以用一個(gè)矩張量表示,即

(3)

式中：M0為等效體力；ni為斷層面法向方向力在x,y,z方向的分量;fi為斷層面滑動方向力在x,y,z方向的分量。因此矩張量又表示為

(4)

如果不考慮標(biāo)量矩M0問題,只考慮點(diǎn)源能量不同方向的變化問題,利用張量研究震源,則為[7]

(5)

如果只考慮3個(gè)主(優(yōu)勢)方向,則為

(6)

本文采用矢量研究微地震震源反演問題。

2 格林函數(shù)正演

2.1 基本方法

矢量反演以地下震源傳播到接收點(diǎn)上的位移能夠被矢量描述為前提,接收點(diǎn)記錄是震源矢量和地層傳播矢量的綜合響應(yīng)結(jié)果。為了進(jìn)行震源反演,必須首先研究正演問題,利用簡化高斯束射線追蹤方法合成地震記錄。微地震監(jiān)測觀測到的微地震記錄是寬角地震記錄,因此我們討論寬角激發(fā)、接收條件下的微地震正演合成記錄。

位移的一般表示式為[8]

(7)

其中

(8)

式中：A(R),A(S)為接收點(diǎn)和源點(diǎn)的振幅;ρ(R),ρ(S)為接收點(diǎn)和源點(diǎn)的密度;v(R),v(S)為接收點(diǎn)和源點(diǎn)的速度;βi,αi是上、下地層P,P′界面兩側(cè)以地層界面的切線方向?yàn)樗捷S的局部坐標(biāo)系中入射線和透射線與局部坐標(biāo)系x正方向的夾角;Ti為透射系數(shù)。

采用Zoeppritz方程計(jì)算透射系數(shù),Zoeppritz方程的矩陣形式為[9]

(9)

式中：α,β,α′,β′分別表示反射P波和反射SV波、透射P波和透射SV波的角度;RPP和RPS是以位移表示的反射P波和反射SV波反射系數(shù);TPP和TPS分別是以位移表示的透射P波和透射SV波透射系數(shù);vP1,vP2,vS1,vS2分別表示下層介質(zhì)的P波速度、上層介質(zhì)的P波速度、下層介質(zhì)的SV波速度、上層介質(zhì)的SV波速度;ρ1,ρ2分別表示上、下層介質(zhì)的密度。

合成地震記錄為

(10)

式中：φ0,φN為中心射線的入射角范圍,包括所有對檢波點(diǎn)R有貢獻(xiàn)的高斯射線束;g(R,φ)是波包;Δφ表示中心射線的角度間隔[10]。

2.2 正演結(jié)果

正演采用均勻水平層狀速度模型,震源子波為雷克子波,10級三分量檢波器接收。為了體現(xiàn)震源能量在觀測位移分量上的變化,設(shè)計(jì)了均勻震源和非均勻震源。正演模擬結(jié)果如圖1所示,在呈現(xiàn)地層響應(yīng)變化的同時(shí),呈現(xiàn)了震源的能量變化引起的位移變化,可見震源非均勻性輻射引起微地震記錄振幅發(fā)生較大變化。

圖1 正演模擬結(jié)果a 均勻點(diǎn)源合成記錄; b 非均勻點(diǎn)源合成記錄

3 微地震矢量反演

3.1 基本方法

反演的基本思想是通過對比觀測的微地震波場運(yùn)動學(xué)和動力學(xué)記錄與理論模型正演波場的運(yùn)動學(xué)和動力學(xué)波場信息,利用迭代方法[11]反演震源參數(shù)。

為了得到準(zhǔn)確的震源矢量反演結(jié)果,建立貼近實(shí)際的初始模型至關(guān)重要。反演模型包括源-檢位置、地層速度模型和震源矢量模型。為了使求解過程快速收斂且得到較高精度的解,提出如下反演策略：①利用走時(shí)反演得到震源位置、計(jì)算格林函數(shù);②對地層速度和震源位置進(jìn)行微調(diào)反演,對震源矢量進(jìn)行調(diào)節(jié)反演。

3.2 基于走時(shí)的距離反演

要獲得震源到接收點(diǎn)的格林函數(shù),需反演得到震源和接收點(diǎn)間的距離[12]。為此,在基于走時(shí)的距離反演基礎(chǔ)上,利用黃金分割法進(jìn)行搜索,并且在對數(shù)域進(jìn)行尋優(yōu)反演,以提高反演精度。具體步驟如下：

1) 設(shè)置計(jì)算區(qū)間[a1,b1]=[Rmin,Rmax],[c1,d1]=[zmin,zmax]及精度要求ε>0,計(jì)算搜索網(wǎng)格點(diǎn)R1=a1+(1-0.618)(Rmax-Rmin),z1=c1+(1-0.618)(zmax-zmin),計(jì)算目標(biāo)函數(shù)值fk,并令k=1。計(jì)算目標(biāo)函數(shù)時(shí)走時(shí)取對數(shù),即fk=lg(|Ti-Tj|)。

2) 如果|fk+1-fk|<ε,則停止計(jì)算,否則轉(zhuǎn)步驟3)。

3) 置ak+1=Rk,ck+1=zk,bk+1=bk,dk+1=dk,ak+1=Rk+1+0.618(bk+1-ak+1),ck+1=zk+1+0.618(dk+1-ck+1),計(jì)算目標(biāo)函數(shù)值fk。如果|fk+1-fk|<ε則停止計(jì)算,否則轉(zhuǎn)步驟4)。

4) 置k=k+1,轉(zhuǎn)步驟2)。

圖2為某一實(shí)際微地震事件的走時(shí)及走時(shí)差,圖3和圖4為利用改進(jìn)后的算法分別在算數(shù)域和對數(shù)域進(jìn)行的距離反演結(jié)果。通過對比可以看出,在徑向方向上,對數(shù)域反演結(jié)果較為收斂和精確。

圖2 各道走時(shí)(a)和相鄰道走時(shí)差(b)

圖3 算數(shù)域反演結(jié)果a 徑向方向; b 縱向方向

圖4 對數(shù)域反演結(jié)果a 徑向方向; b 縱向方向

3.3 震源矢量反演

3.3.1 震源矢量模型建立

利用偏振分析技術(shù),通過觀測的三分量記錄得到特征向量進(jìn)而得到總的觀測主矢量;利用格林函數(shù)正演得到地層速度、源-檢位置格林函數(shù)矢量。將偏振分析得到的矢量去除格林矢量后作為震源初始矢量模型。

3.3.2 反演尋優(yōu)方法

利用模擬退火方法[13]進(jìn)行反演,參與反演調(diào)節(jié)的模型為矢量模型。

3.3.3 目標(biāo)函數(shù)的確定

三分量檢波器坐標(biāo)雖然利用射孔資料進(jìn)行了標(biāo)定,但是由于射孔源的非均勻性、地層效應(yīng)、檢波器響應(yīng)及其它干擾因素的影響,標(biāo)定結(jié)果很不準(zhǔn)確。所以,作為三分量振幅反演,不但要進(jìn)行速度模型調(diào)節(jié)、震源輻射能量[14]調(diào)節(jié)反演,還要進(jìn)行三分量檢波器坐標(biāo)的標(biāo)定調(diào)節(jié)反演。只有在檢波器坐標(biāo)確定的條件下,才能得到振幅在其上的準(zhǔn)確投影分量。格林函數(shù)地層響應(yīng)的極化矢量為AG,通過格林函數(shù)正演求得,震源的極化矢量為AS,總的極化矢量為A=AG+AS,則理論模型計(jì)算的3個(gè)分量為

(11)

式中：φ0是極化矢量的傾角,是確定的,因?yàn)槿至繖z波器垂直分量的方向是固定的;α0是檢波器水平分量的初始方位角,與射孔標(biāo)定精度有關(guān),是不確定的。α=αs-α0,αs是事件s的方位角。

由于同一個(gè)事件對于多級檢波器的絕對方位角(事件方位角與檢波器初始方位角的差)是相同的,因此,反演目標(biāo)函數(shù)加入了多級檢波器絕對方位角相等的約束條件,以使反演解得到更好的控制。即

(12)

式中：i是檢波器級數(shù),j是相鄰事件個(gè)數(shù)序號。

最終反演目標(biāo)函數(shù)為

(13a)

αij=αs,i,j-α0,j=αi+1,j=αs,i+1,j-α0,j

i=1,2,…,N-1

j=1,2,…,K

(13b)

式中：Ux,g,Uy,g,Uz,g表示實(shí)際觀測的3個(gè)分量；p為模數(shù)。

為了建立三分量檢波器較準(zhǔn)確的坐標(biāo)系,利用多個(gè)射孔資料進(jìn)行矢量反演,具體步驟如下：

1) 利用走時(shí)距離反演結(jié)果及建立的速度模型,進(jìn)行格林函數(shù)正演;

2) 令k=1;

3) 將原走時(shí)距離反演建立的速度模型作為初始速度模型,源矢量反演時(shí)對其進(jìn)行微調(diào),即v=v0+Δv;

4) 震源矢量調(diào)節(jié),對觀測的三分量數(shù)據(jù)進(jìn)行偏振分析,將求得的偏振矢量作為矢量反演的初始矢量模型,在此基礎(chǔ)上進(jìn)行擾動調(diào)節(jié),即As=As0+ΔAs;

5) 檢波器初始方位角調(diào)節(jié),將射孔資料計(jì)算的檢波器方位角作為初始方位角,在反演過程中進(jìn)行不斷調(diào)節(jié),即α0=α0,s+Δα0;

7)k=k+1,轉(zhuǎn)向步驟3);

8) 停止計(jì)算。

4 理論模型驗(yàn)證分析

采用均勻水平層狀模型驗(yàn)證了震源矢量反演方法的正確性。設(shè)置模型事件方位角為60°,震源位置(200m,-200m,1000m),井下垂直10級三分量檢波器接收,第1級檢波器接收位置(0,0,950m),各級檢波器之間的距離為20m。為了驗(yàn)證源矢量方向和方位方向的不一致,采用能量大小不同、方向不同的3個(gè)源矢量合成一個(gè)總的源矢量,合成矢量方向方位角見表1。對合成記錄進(jìn)行偏振分析求得各級檢波器方位角,再進(jìn)行震源矢量反演,結(jié)果如表1所示。分析模型和反演結(jié)果可見,通過合成記錄計(jì)算的方位角和模型方位角差距較大,震源反演的方位角基本上接近模型合成矢量方位角。說明震源矢量反演方法正確,反演結(jié)果具有較高的精度。

表1 均勻水平層狀模型分析結(jié)果 °

5 實(shí)際資料驗(yàn)證分析

利用某油田壓裂微地震監(jiān)測資料進(jìn)一步驗(yàn)證震源矢量反演方法的正確性。該區(qū)壓裂井和監(jiān)測井兩井之間的距離為657.8m,目的層段深度2130.0m,壓裂巖石巖性為灰?guī)r。采用12級三分量檢波器接收,各級檢波器之間的距離為10m,硬鏈接3級,間距為4.5m。圖5a為直接反演定位結(jié)果,圖5b為震源矢量反演方法定位結(jié)果。

研究區(qū)壓裂目的層巖性較均勻,壓裂裂縫的發(fā)育主要受控于水壓力和地應(yīng)力圍壓分布,主構(gòu)造發(fā)育方向近似為東北方向;又由井資料可知,目的層最大主應(yīng)力方向?yàn)闁|北方向。據(jù)此推斷,在水力壓力和地應(yīng)力作用下,巖石破裂裂縫的延伸方向應(yīng)該為東北方向。該區(qū)壓裂射孔點(diǎn)平面位置在(257m,10m)處,對比分析圖5a和圖5b可見,直接反演定位結(jié)果的射孔點(diǎn)附近事件分布稀疏,而射孔點(diǎn)附近以外事件分布較稠密,不符合壓裂微地震事件發(fā)育分布規(guī)律;通過震源矢量反演再定位后,射孔點(diǎn)附近事件分布稠密,而射孔點(diǎn)附近以外事件分布稀疏,這與壓裂微地震事件發(fā)育分布規(guī)律相符。

圖5 某油田壓裂微地震定位結(jié)果a 直接反演定位結(jié)果; b 震源矢量反演定位結(jié)果

6 結(jié)束語

微地震震源矢量反演方法定位結(jié)果較直接反演方法定位結(jié)果在精度和可靠程度方面有了較大提高。這是因?yàn)椋孩倮煤喕母咚故渚€追蹤方法合成地震記錄,走時(shí)和振幅精度能夠滿足反演的要求;②改進(jìn)了走時(shí)距離反演,將以往的常規(guī)網(wǎng)格搜索方法改進(jìn)為黃金分割網(wǎng)格搜索方法;③在對數(shù)域進(jìn)行走時(shí)正演和反演迭代計(jì)算,提高了徑向距離對走時(shí)的敏感性;④在進(jìn)行速度模型反演調(diào)節(jié)、震源輻射能量反演調(diào)節(jié)時(shí),利用多個(gè)射孔資料進(jìn)行三分量檢波器坐標(biāo)的標(biāo)定調(diào)節(jié),得到了相對準(zhǔn)確的檢波器坐標(biāo)方位;⑤根據(jù)同一個(gè)事件對于多級檢波器的絕對方位角相同,反演目標(biāo)函數(shù)增加了多級檢波器絕對方位角相等的約束條件,使反演結(jié)果更加可靠。

參 考 文 獻(xiàn)

[1] 宋維琪,陳澤東,毛中華.水力壓裂裂縫微地震監(jiān)測技術(shù)[M].北京：中國石油大學(xué)出版社,2008:156-167

Song W Q,Chen Z D,Mao Z H.Hydro-fracturing break microseismic monitoring technology[M].Beijing:China University of Petroleum Press,2008:156-167

[2] 許沖,徐錫偉,于貴華.基于證據(jù)權(quán)方法的玉樹地震滑坡危險(xiǎn)性評價(jià)[J].地震地質(zhì),2013,35(1):151-164

Xu C,Xu X W,Yu G H.The Yushu earthquake triggered landslide hazard evaluation based on weight of evidence method[J].Seismology and Geology,2013,35(1):151-164

[3] 謝飛,李佩,黃中玉,等.高斯射線束疊前深度偏移成像研究[J].石油物探,2013,52(1):65-71

Xie F,Li P,Huang Z Y,et al.Gaussian beam prestack depth migration[J].Geophysical Prospecting for Petroleum,2013,52(1):65-71

[4] 黃建平,張晴,張凱,等.格林函數(shù)高斯束逆時(shí)偏移[J].石油地球物理勘探,2014,49(1):101-106

Huang J P,Zhang Q,Zhang K,et al.Reverse time migration with Gaussian beams based on the Green function[J].Oil Geophysical Prospecting,2014,49(1):101-106

[5] Engell-Sorensen L.Inversion of arrival times of microearthquake sources in the North Sea using a 3-D velocity structure and prior information,part I:method[J].Bulletin of the Seismological Society of America,1991,81(4):1183-1194

[6] Oladapo M I.Linearization of Zoeppritz equations and practical utilization[J].International Journal of Physical Sciences,2013,8(24):1298-1306

[7] Prugger A F,Gendzwill D J.Microearthquake loca-tion:a nonlinear approach that makes use of a simplex stepping procedure[J].Bulletin of the Seismological Society of America,1988,78(2):799-815

[8] 安藝敬一,理查德.定量地震學(xué):理論與方法[M].北京:地震出版社,1980:365-398

Aki K,Richards P G.Quantitative seismology:theory and methods[M].Beijing:Seismological Press,1980:365-398

[9] McMechan G A.Migration by extrapolation of time-dependent boundary values [J].Geophysical Prospecting,1983,31(3):413-420

[10] Furumura T,Kennett B,Takenaka H.Parallel 3-D pseudospectral simulation of seismic wave propagation[J].Geophysics,1998,63(1):279-288

[11] Blake B,Figueroa D,Hofland G,et al.3D forward ray trace seismic modeling of strike lines in complex geology[J].Expanded Abstracts of 69thAnnual Internat SEG Mtg,1999,1871-1874

[12] Kummer B,Behle A,Dorau F.Hybrid modeling of elastic-wave propagation in two-dimensional laterally inhomogeneous media[J].Geophysics,1987,52(6)：765-771

[13] Stoffa P L,Sen M K.Nonlinear multiparameter optimization using genetic algorithms:inversion of plane-wave seismograms[J].Geophysics,1991,56(11):1794-1810

[14] Pei D,Quirein J A,Cornish 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