宋維琪，朱海偉，姜宇東，郭全仕，曹 輝

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

微地震監(jiān)測研究的重點(diǎn)在于微地震事件的準(zhǔn)確定位。影響精確定位的因素很多，主要包括速度模型建立精度、反演算法的適用性、正演算法的精度等方面［1－3］。高精度的反演研究，目前還主要以時差反演為主，通過模型校正和反演方法的研究來提高定位精度，其中又分為P波時差反演和P－S波時差反演［4－5］。

常規(guī)地震反演中，共軛梯度法、模擬淬火法、遺傳算法等線性或非線性算法都能夠取得較好的反演效果。但是，微地震反演中，由于信號的信噪比低，加上淺地表各種噪聲的干擾，難以拾取準(zhǔn)確的初至，使得單純地移植這些反演方法無法取得滿意的效果。一般來講，線性反演方法速度快，效率高，但是易于局部收斂，對設(shè)定的初值要求高。非線性方法全局搜索能力強(qiáng)，但是算法設(shè)計比較復(fù)雜，反演參數(shù)越多，反演速度越慢［6］。

對地面微地震反演，由于實(shí)際接收到的檢波器信號中有些道的同相軸信息基本缺失，使得能夠用于反演的道非常少；其次，對于變量的一階導(dǎo)數(shù)的求取也非常困難，主要是由于目標(biāo)函數(shù)值對于參量的擾動非常不敏感，使得求導(dǎo)運(yùn)算十分復(fù)雜。此外，不同變量的擾動，對于目標(biāo)函數(shù)值的影響也不同，例如對于相同的變化量Δx，變量層速度v或水平距離r對目標(biāo)函數(shù)的影響完全不同，使得迭代過程十分繁瑣。

遺傳算法是一種非線性反演方法，其獨(dú)有的個體、種群概念保證了搜索的路徑是多維的，可以避免單維搜索陷入局部極值。但是，遺傳操作過程中主要算法參數(shù)的選擇對于反演的影響很大，既影響反演速度，又影響反演精度。也就是說，在微地震事件反演中，如果任意定義反演的搜索范圍，任意定義選擇算子、交叉算子、變異算子、適應(yīng)度函數(shù)、種群大小、遺傳代數(shù)等，那么遺傳算法反演的效率與精度都會發(fā)生很大的變化，也增大了反演結(jié)果的不穩(wěn)定性。

地面微地震監(jiān)測星形的觀測方式和資料信噪比低的特殊性［7－8］，尤其是由于有用信號十分微弱，使得拾取的初至信息不準(zhǔn)確，以往確定性的線性或非線性反演方法幾乎不能取得可靠的反演結(jié)果。貝葉斯反演方法是把先驗(yàn)信息和后驗(yàn)信息結(jié)合起來的推理估計方法，可以更好地融合多方面先驗(yàn)信息實(shí)現(xiàn)后驗(yàn)信息的最佳估計。因此，結(jié)合地面微地震資料的特點(diǎn)，研究了地面微地震資料震源定位的貝葉斯反演方法。

1 方法原理

地面微地震監(jiān)測采用多線觀測，利用多線觀測結(jié)果進(jìn)行微地震事件的定位。根據(jù)全概率理論，可以把最終定位解做為事件A，每條線都以不同概率確定事件A。把各條線反演問題作為最終反演結(jié)果的一個劃分即Bi（i＝1，2，…，M），其全概率公式為

式中：P（Bi）為先驗(yàn)概率，根據(jù)每條線的信噪比確定；P（A｜Bi）為后驗(yàn)概率。

貝葉斯反演理論把反演問題與觀測數(shù)據(jù)信息、模型信息以及先驗(yàn)信息通過貝葉斯推理聯(lián)系起來。貝葉斯反演有多種求解方法，其中最常用的是貝葉斯最大后驗(yàn)解方法。設(shè)G為微地震到時觀測數(shù)據(jù)，X為待反演參數(shù)，模型誤差與待反演參數(shù)的先驗(yàn)估計均為正態(tài)分布，對每條線反演求解，建立目標(biāo)函數(shù)求取貝葉斯最大后驗(yàn)解［9］：

式中：F（X）為模型結(jié)果；XP為參數(shù)的先驗(yàn)估計值；CD和CM為模型誤差與先驗(yàn)估計誤差的協(xié)方差矩陣；λ為權(quán)重系數(shù)，λ＝CM／CD。目標(biāo)函數(shù)S（X）的求解，就是在數(shù)據(jù)擬合（第1項(xiàng)）和先驗(yàn)估計（第2項(xiàng)）之間尋求一種折衷。權(quán)重系數(shù)越大，反演解越接近先驗(yàn)估計；反之，反演解越接近于不穩(wěn)定的原問題。非線性反演問題權(quán)重系數(shù)的選取是一個廣受關(guān)注的問題，在無法確定觀測數(shù)據(jù)噪聲的情況下，比較常用的有廣義交叉驗(yàn)證（Generalized Cross－Validation，GCV）方法和L－曲線（L－Curve）方法等。我們采用赤池貝葉斯信息準(zhǔn)則（Akaike's Bayesian information criterion，ABIC）來確定最佳權(quán)重系數(shù)。ABIC是基于相對信息熵理論而提出的權(quán)重系數(shù)選取方法，其最佳權(quán)重系數(shù)由最小化ABIC獲得。ABIC的表達(dá)式為

式中：N為CD的維數(shù)；J為雅可比矩陣；c為常數(shù)。實(shí)際求解時，首先選定若干λ值并求得相應(yīng)的反演解，然后根據(jù)λ與反演解的值計算ABIC，ABIC 最小值對應(yīng)的λ即為最佳權(quán)重系數(shù)。假定抽樣計算獨(dú)立，則評價函數(shù)如

式中：σD，σM分別為理論模型和觀測值殘差的方差、模型參數(shù)的方差。這時，λ＝σD／σM。

對參數(shù)（速度、震源位置）的最優(yōu)估計問題，等價為在貝葉斯理論下，具有先驗(yàn)信息條件時，經(jīng)過貝葉斯推理，得到目標(biāo)函數(shù)的極大后驗(yàn)概率。從公式中看到，目標(biāo)函數(shù)的極大后驗(yàn)概率是理論模型和觀測值殘差、模型參數(shù)的參差極大后驗(yàn)概率的加權(quán)組合。實(shí)現(xiàn)時，通過抽樣計算得到參數(shù)與速度參數(shù)模擬參差結(jié)果的先驗(yàn)分布p（Y，θ），確定θ的后驗(yàn)密度ρ（θ，Y）。設(shè)h（θ，Y）為θ和Y的聯(lián)合概率密度函數(shù)，則

式中：m（Y）為Y的邊緣密度函數(shù)，定義為

由（3）式和（4）式得

式中：p（Y）為樣本Y＝（y1，y2，…，yn）的聯(lián)合概率密度函數(shù)。

可見，θ的后驗(yàn)密度ρ（θ，Y）估計取決于先驗(yàn)分布ρ（θ）和已知θ情況下的樣本的條件概率密度函數(shù)p（Y，θ）（或稱為已知Y條件下，θ的似然函數(shù)）。p（Y，θ）通常采用模型的計算值與實(shí)測值殘差的函數(shù)或者正態(tài)分布表示［10－11］，如

式中：σ2為殘差的方差。

理論正態(tài)分布似然函數(shù)表示為

后驗(yàn)概率表示為

分析（10）式可知，在給定樣本空間情況下，決定極大后驗(yàn)概率的是概率密度函數(shù)和分布函數(shù)的積分值。如果數(shù)據(jù)出現(xiàn)奇異，則積分結(jié)果奇異。因此，為了獲得穩(wěn)定的目標(biāo)函數(shù)值，需要從計算的殘差方差中利用回歸方法或擬合方法得到理論模擬概率密度函數(shù)，然后計算后驗(yàn)概率，最后從不同速度模型對應(yīng)的后驗(yàn)概率中取極大后驗(yàn)概率為目標(biāo)函數(shù)結(jié)果。

2 尋優(yōu)方法

公式（1）中右邊第1項(xiàng)是一個有關(guān)地層速度和震源位置的多模態(tài)函數(shù)［12－13］，需要采用非線性法全局搜索算法來優(yōu)化。為提高計算效率，采用極快速模擬退火方法加網(wǎng)格法的混合算法作為搜索方法。首先利用網(wǎng)格算法使搜索落入最優(yōu)解所在的凸區(qū)間，再采用極快速模擬退火算法搜索最優(yōu)解。這樣既可以防止算法收斂于局部極值點(diǎn)，又極大地提高了算法的收斂速度。

3 建模方法

地面微地震走時記錄由中深層速度變化較緩部分和淺地表速度變化相對劇烈部分組成。在反演時這兩部分應(yīng)分別對待，淺部采用橫向變速模型，中深部采用水平橫向均勻速度模型。在反演過程中，中深部數(shù)據(jù)和淺部數(shù)據(jù)可能會相互抑制，即淺部數(shù)據(jù)如果擾動劇烈，可能會把中深層的數(shù)據(jù)緩慢擾動淹沒，使中深部反演變得遲鈍。同樣，中深部反演數(shù)據(jù)的擾動變化平緩，會影響淺部反演的敏感性。為此，反演迭代過程也分為兩部分，即

式中：tl為理論計算的走時；tq為理論計算的淺部走時；ts為理論計算的深部走時，通過射線追蹤正演計算。

計算淺部理論到時采用的近似方法為：假設(shè)淺部地層速度較低，地震波的射線路徑與水平地面垂直，不考慮橫向速度變化引起的射線路徑的彎曲變化。

式中：hj為測井起始深度到高程參考基準(zhǔn)面的厚度；hgc（x）為高程；vq為淺部速度模型。

目標(biāo)函數(shù)（2）式相應(yīng)地變化為

4 利用井資料和射孔資料反演初始速度模型

微地震壓裂監(jiān)測過程中，通常事先都要在壓裂井段進(jìn)行射孔以獲取射孔事件記錄。通過射孔記錄到時和射孔點(diǎn)的已知位置進(jìn)行速度模型的校正和建立，進(jìn)而獲取更加準(zhǔn)確的初始速度模型。地面微地震監(jiān)測不同于井中微地震監(jiān)測，一般沒有淺部測井資料。因此，地面微地震射孔資料不單是用來校正速度模型，還要用來進(jìn)行淺部速度模型的建立。中深層初始速度模型資料一般由測井資料獲得，而淺部速度資料通常難以收集，因此通過射孔反演方法求取淺部速度。求取步驟如下：

1）利用井資料建立中深層速度模型；

2）計算中深層模型的走時；

3）計算淺部走時

4）求淺部初始速度。

假定微地震波走時路徑垂直水平面，則各接收點(diǎn)到參考面的初始速度為

至此，就建立了從淺部到深部整個速度初始速度模型。然后采用本文研究的貝葉斯反演方法進(jìn)行速度模型的反演。具體的反演步驟：

1）在初始速度模型上進(jìn)行速度參數(shù)的擾動調(diào)節(jié)

并計算其方差。

2）進(jìn)行正演計算理論走時，然后計算觀測到時和理論計算到時的殘差與方差以及加權(quán)系數(shù)λq，λs。

3）對速度進(jìn)行參數(shù)隨機(jī)抽樣，擬合（或回歸分析）出一個理論殘差方差的后驗(yàn)概率密速函數(shù)，即目標(biāo)函數(shù)的后驗(yàn)概率密度函數(shù)、加權(quán)函數(shù)后驗(yàn)密度函數(shù)、速度參數(shù)方差的后驗(yàn)概率密度函數(shù)。通過以上各變量的函數(shù)求得雅可比矩陣。根據(jù)（15）式計算評價函數(shù)的極小值。此評價函數(shù)的極小值對應(yīng)的速度參數(shù)就是反演結(jié)果。

4）反演每一條測線，得到了各線的淺地表速度模型和中深部地層速度模型。

5 微地震事件震源位置反演

通過測井資料和射孔資料反演出各條線整個地層的速度模型后，根據(jù)拾取的走時反演微地震震源位置。根據(jù)反演要求的不同，可以分為位置參數(shù)反演和速度、位置參數(shù)同時反演，我們只考慮位置參數(shù)反演。利用拾取的各條線到時和上面建立的速度模型，采用本文反演方法進(jìn)行反演，得到震源位置。由于各條測線微地震記錄的信噪比不同，拾取的到時和建立的速度模型精度不同，反演結(jié)果也不盡相同。綜合各條線的反演結(jié)果，根據(jù)全概率公式推斷估計最佳震源位置。

5.1 模型試驗(yàn)

設(shè)單個理論模型震源坐標(biāo)為（200m，300m，2 000m），采用上文建立的速度模型進(jìn)行反演。理論模型的走時利用射線追蹤正演獲得，并在理論走時中加上一定的擾動噪聲。圖1為我們研究的貝葉斯方法反演結(jié)果和以前方法反演結(jié)果的對比，可以看出，貝葉斯反演方法具有較好的抗噪能力。

圖1 理論模型資料用以前方法反演結(jié)果（a）和貝葉斯方法反演結(jié)果（b）對比

5.2 實(shí)際資料反演

對實(shí)際地面微地震資料進(jìn)行一定處理后，拾取有效事件的初至，結(jié)果如圖2 所示。然后對利用射孔資料反演校正的速度模型進(jìn)行反演。反演結(jié)果如圖3所示。圖3中紅色和紫色的實(shí)點(diǎn)表示反演的2 個事件位置，淺藍(lán)色的圈表示射孔點(diǎn)。反演結(jié)果和國外商業(yè)軟件反演結(jié)果基本一致，如表1所示。

表1 2個實(shí)際微地震事件位置反演結(jié)果比較

6 結(jié)束語

地面微地震震源反演問題由于各種因素的影響和資料信噪比低，具有較大程度的不確定性。因此，將先驗(yàn)信息和后驗(yàn)信息結(jié)合起來進(jìn)行最優(yōu)解的評價，較好地實(shí)現(xiàn)了不確定性反演問題的尋優(yōu)過程。將反演建模和迭代過程分為兩部分進(jìn)行，較好地解決了淺部和深部速度變化程度不同對反演解的影響。采用極快速模擬退火方法加網(wǎng)格法的混合算法作為搜索方法，網(wǎng)格法可以使搜索過程落入最優(yōu)解所在的凸區(qū)間，而極快速模擬退火算法搜索最優(yōu)解，既可以防止算法收斂于局部極值點(diǎn)，又能極大地提高算法的收斂速度。反演結(jié)果表明，反演方法在抗擾動噪聲和解決不確定性方面具有顯著優(yōu)勢?；谧畲蠛篁?yàn)估計的貝葉斯反演方法對先驗(yàn)信息估計模型具有非常大的依賴性，因此，下一步對先驗(yàn)信息估計模型建立還需進(jìn)行深入研究。

［1］Stoffa P L，Sen M K.Nonlinear multiparameter optimization using genetic algorithms：inversion of planewave seismograms［J］.Geophysics，1991，56（11）：1794－1810

［2］Pei D H，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）：47－55

［3］Lisbeth E S.Inversion of Arreal times of microearthquake sources in the noth sea using a 3－D velocity structure and prio information：part I［J］.Method Bulletin of the Seismological Society of America，1991，81（4）：1183－1194

［4］Willis R B，F(xiàn)ontaine J，Paugh L，et al.Geology and geometry：a review of factors affecting the effectiveness of hydraulic fracturing［R］.SPE Eastern Regional Meeting，97993－MS，14－16

［5］Chunduru R K，Sen M K，Stoffa P L.Hybrid optimization methods for geophysical inversion［J］.Geophysics，1997，62（4）：1196－1207

［6］Phillip W S，House L S.Micro－seismic mapping of a Cotton Valley Hydraulic Fracture using decimated downhole arrays［J］.International Exposition and Sixty Eighth Annual Meeting，1998，9：13－18

［7］Dumay J，F(xiàn)ournier F.Multivariate statistical analyses applied to seismic facies recognition［J］.Geophysics，1988，53（9）：1151－1159

［8］De Meersman K，Kendall J M，van der Baan M.The 1998 Valhall microseismic data set：an integrated study of relocated sources，seismic multiplets，and Swave splitting［J］.Geophysics，2009，74（6）：183－195

［9］詹海剛，施平，陳楚群.基于貝葉斯反演理論的海水固有光學(xué)特性準(zhǔn)分析算法［J］.科學(xué)通報，2006，51（2）：204－210 Zhan H G，Shi P，Chen C Q.Based on Bayesian inversion theory of inherent optical properties of seawater anaysis algorithm［J］.Chinese Science Bulletion，2006，512（2）：204－210

［10］王華忠，方正茂，徐兆濤，等.地震波旅行時計算［J］.石油地球物理勘探，1999，34（2）：155－163 Wang H Z，F(xiàn)ang Z M，Xu Z T，et al.Computation of seismic travel time［J］.Oil Geophysical Prospecting，1999，34（2）：155－163

［11］萬永革，李鴻吉.遺傳算法在確定震源位置中的應(yīng)用［J］.地震地磁觀測與研究，1995，16（6）：1－7 Wan Y G，Li H J.The preliminary study on the seismic hypocenter location using genetic algorithms［J］.Seismological and Geomagnetic Observation and Research，1995，16（6）：1－7

［12］屈永華，王錦地，劉素紅，等.貝葉斯網(wǎng)絡(luò)支持的地表參數(shù)混合反演模式研究［J］.遙感學(xué)報，2006，10（1）：6－13 Qu Y H，Wang J D，Liu S H，et al.Study on hybrid inversion scheme under Bayesian network［J］.Journal of Remote Sensing，2006，10（1）：6－13

［13］宋維琪，劉軍，陳偉.改進(jìn)射線追蹤算法的微震源反演［J］.物探與化探，2008，32（3）：275－277 Song W Q，Liu J，Chen W.Microearthquake source inversion of an improved ray tracing algorithm［J］.Geophysical and Geochemical Exploration，2008，32（3）：275－277