李自红,陈 慧

(1.山西省地震局,太原 030021;2.太原理工大学 地震与地质灾害防治研究所,太原 030024;3.太原大陆裂谷动力学国家野外科学观测研究站,太原 030025)

我国是世界上发生煤矿灾害事故最严重的国家之一,随着煤矿开采深度的增加和开采强度的提高,煤与瓦斯突出、顶板垮落、冲击地压等煤岩动力灾害日趋严重[1]。冲击地压是采煤诱发的破坏性矿震,是威胁煤矿安全、高效生产的重大灾害之一。而煤层气(瓦斯)是在成煤作用过程中相伴而生成,是一种自生自储式非常规天然气。实践表明,地面煤层气井采前预抽是低透气性煤层瓦斯灾害治理最有效的方式之一,而水力压裂是当前煤层气井增产的有效措施[2]。

无论是采煤引发的冲击地压、还是煤层气开采过程的水力压裂,都伴随着众多小能量级别的微震,而这些微震活动则是研究冲击地压或评价水力压裂效果的重要信息。通过微震监测,可以记录煤岩体破裂过程中产生的微震信号,并对微震震源进行精确定位。定位结果不但可以用于预测冲击地压发生的位置、时间和规模,也是评价水力压裂后的煤层裂缝形态和展布规律的主要参数。

震源位置是微震监测中需要确定的最关键和最基本的参数之一[3]。高精度的震源位置确定不仅依赖于良好的震相数据,合理的速度模型和可靠的定位方法也是必不可少的。国内外众多学者一直把提高地震定位精度作为一个重要的研究课题,并不断改进或提出新的地震定位方法[4]。早期的地震定位方法(几何作图法)主要是以走时方程为理论依据,根据台站位置及走时资料,通过作图确定震源位置,其精度较低[5]。随着计算机技术在地震定位中的广泛应用,不仅常规地震定位方法得到不断改进,适合冲击地压、水力压裂监测等微震定位方法也得到迅速发展。

1 常规地震定位方法

地震定位的实质是根据地震事件到各个台站的震相到时来反演地震的震源参数(纬度、经度、深度和发震时刻)。常用的地震定位方法有几何作图法、线性定位法(包括绝对定位法和相对定位法)及非线性定位法。

1.1 绝对定位法

上世纪初德国物理学家Geiger[6]提出一种单地震事件绝对定位方法,是经典定位方法。其原理是将含震源时空参数的非线性方程作一阶泰勒展开,将非线性走时方程线性化,并通过最小二乘法求解参数方程来获取震源参数。20世纪70年代随着计算机技术的发展,Geiger的经典定位方法才得以应用。后人以经典定位法为基础,提出了一系列的改进的经典定位法。Lee等[7]在J.Eaton的HYPOLAYR[8]程序基础上于上世纪七十年代初期完成了HYPO71地震定位程序,此后相继给出了HYPO78—86系列程序,现已形成HYPOELLIPSE定位程序[9]。Klein[10]提出HYPOINVERSE,Lienert等[11]在此基础上提出可以定位远震的HYPOCENTER,Nelson和Vidale[12]也改进了HYPOINVERSE,提出了三维速度模型下的QUAKE3D方法。

Geiger经典定位方法对初始条件的依赖性较大,有时因台站分布不合理,地震偏在台网的一侧等,会使得定位发散或得不到定位结果[13]。为了解决这些问题,人们又提出了一些改进方法,如迭代收敛法、阻尼最小二乘法、共轭梯度法以及奇异值分解法等,在一定程度上解决了迭代过程的失稳和发散,显著提高了数值计算的稳定性[14]。

然而经典定位法由于受复杂的地壳结构影响,其定位精度有限。Douglas[15]在1967年提出了震源位置与台站校正联合反演的JED方法。后来Dewey、Pavlis及Pujol等人相继对此方法做出改进、简化,使其变为现在应用较为广泛的JHD定位方法[16-18]。王椿镛等[19]根据昆明台网区域地震资料,用JHD和参数分离法对各台站P波走时进行校正,使定位精度有了较大提高。

Crosson[20]提出震源位置和速度结构的联合反演理论(即SSH方法),可以削弱地震定位中速度模型引起的误差。以此为基础,一系列的震源位置与速度结构反演方法被提出,如Aki等[21]提出的三维速度结构与震源位置联合反演方法等。Kissling[22]等通过研究震源位置和一维速度模型耦合问题,提出了确定最佳一维速度模型的Velest方法。该方法计算的速度模型可使定位结果的走时残差均方根最小,定位精度得到显著提高,被国内外学者广泛地应用于地震精度定位初始速度模型的确定。

1.2 相对定位法

为了减小速度模型对定位结果的影响,不少学者提出了使用到时差的相对定位法,包括主事件定位法(ATD)和双差定位法(HYPODD),是单纯使用震相数据就可以获得高精度地震参数的方法。

Spence[23]在JED定位方法的基础上提出主事件定位法(ATD)。其基本原理是选定一个震源位置比较精确的地震作为主事件,通过计算其周围一群地震事件相对于它的位置来确定震源位置,缺点是不能应用于空间跨度较大的地震事件群体定位,另外也不适合近震定位。周仕勇等[24]对主事件定位方法作了较大改动,定位方程中没有包含待定地震发震时刻的未知量,发震时刻未知量是在确定待定地震震源位置后,根据波的传播距离和传播速度计算,避开了发震时刻的直接求解,并且釆用首波到时资料对震源深度专门进行确定。

Waldhauser和Ellsworth[25]提出了双差定位法(HYPODD)。其基本原理是,对于发生在同一区域的多个地震事件的地震序列,选取满足一定条件的两个地震事件组成地震对,利用地震对的观测走时差减去理论计算值的走时差的残差(简称“双差”)确定每个地震事件的相对位置,并在反演中引入“质心”,表示所有的地震经重新定位后其平均“位移”为0,以约束反演结果。Zhang、Thurber[26]在双差定位法的基础上提出了双差层析成像方法,该方法使用水平层状速度模型对震群重定位的同时能够得到改进的速度结构。

1.3 非线性定位方法

地震定位的实质是根据一个或多个地震事件到各个台站的震相到时来反演地震的震源参数,本来是一个非线性问题。除前面提到的线性定位方法外,也可以在震源参数空间内使用最优化方法直接求解,即所谓的非线性定位方法。目前常用的非线性定位法有Powell法和遗传算法,另外还包括模拟退火法、单纯形法、蒙特卡罗法等。非线性定位方法因运算量大、效率较低,目前使用较少。随着计算机性能的提高,非线性定位方法将可能会成为普遍使用的常规地震定位方法。

2 无需预先测速的微震震源定位方法

纵观国内外广泛使用的地震及微震定位方法,都以预先给出速度模型为前提,速度模型的准确性直接影响着定位精度。为解决速度模型影响定位精度的问题,基于X. B. Li和L. J. Dong[27]无需测速的二维定位方法,董陇军给出了无需预先测速的微震三维定位方法的数学拟合形式,按其因变量为到时、到时差、到时差商分别给出了计算公式[28]。

2.1 因变量为到时的拟合形式

假定微震震源到各台站间为均匀速度模型且传播速度未知,将其用c表示;震源坐标为(x0,y0,z0);Ti(i=1,2,L,…,n)为第i个台站,各台站坐标为(xi,yi,zi)(i=1,2,L,…,n);li(i=1,2,L,…,n)为各台站至震源的距离;t0为微震发生的时刻,ti(i=1,2,L,…,n)为各台站记录的P波到达时刻。则有

(1)

台站至震源的距离li为

(2)

将式(2)代入式(1)中,可得第i个台站计算到时为

(3)

式中:(xi,yi,zi),ti为已知量;P波的传播速度c;微震的震源位置(x0,y0,z0)和发生时刻t0为未知量。

设l为各台站至震源的距离平均值,t为P波到达各台站的时刻平均值,则有

(4)

其中,

(5)

由式(3),(4)可以构成最小二乘函数:

(6)

式(6)是一个非线性拟合问题,对其进行最小二乘求解,便可得到震源位置(x0,y0,z0)、微震发生时刻t0及介质传播速度c的解。该方法通过4个已知参量拟合5个未知参量,因此至少需要布设5个台站。

2.2 因变量为到时差的拟合形式

假定微震震源到各台站间为均匀速度模型且传播速度未知,将其用c表示,第k个台站计算到时可表示为

(7)

2个台站i和j的到时之差为

(8)

对于每一组观测值(xik,yik,zik;xjk,yjk,zjk),式(8)可确定一个拟合值:

(9)

因式(10)为x0,y0,z0,c的二次非负函数,因此其最小值总是存在的,故总可以得到震源坐标(x0,y0,z0)与介质传播速度c。如只解决震源定位问题,只需拟合式(10)的x0,y0,z0即可。该方法通过4个已知参量拟合4个未知量,因此至少需要布设4台站。

2.3 因变量为到时差商的拟合形式

对于因变量为到时差商的拟合形式,其拟合量只有震源位置(x0,y0,z0)。但因其分母项为两个台站到震源距离之差,如果出现两个或以上台站到震源距离相等的情况,这些台站将不会参与震源定位,因此会影响到定位精度。为提高定位精度,该方法比上述其他两种方法需要布设更多的台站,成本较为昂贵。其数学形式可参考董陇军等有关方面的文献。

3 改进的震源扫描定位方法

3.1 震源扫描算法的原理

本世纪初Kao和Shan提出震源扫描算法(Source-ScanningAlgorithm,SSA)[29],是一种对震源分布进行成像的方法,该方法充分利用地震波形资料,在不用精确拾取震相到时和计算理论走时的情况下达到比较理想的定位效果[30]。通过在整个可能的时间、空间内搜寻,给出最可能的震源位置,其基本原理如下。

假设N个台站记录到某一微震事件,首先对每个台站记录的地震波形进行归一化处理,然后计算某个点η在某个时刻τ的“亮度”,其计算公式为

(11)

式中:N为台站个数;μn为第n个台站记录到的地震波形归一化结果;tη n为从点η处到台站n的走时;τ是微震事件的发生时刻。如果所有的最大振幅都是由点η和时间τ产生,则br(η,τ)=1,但br(η,τ)=1的情况很少见。可以认为整个时空中“亮度”函数取最大时的点η和时间τ,分别是微震的震源位置和发生时刻。

上述问题需要利用网格搜索算法求解,通过遍历所有网格点来寻找整个空间中的最大“亮度”值。其精度主要依赖于网格的尺度大小,尺度越大反演精度越低,尺度越小反演精度越高、但计算量也越大。微震定位中待优化的参数包括介质传播速度V、震源位置(x0,y0,z0)及发震时刻t0等;进行网格搜索时,每个参数都需给出解的范围,并按一定尺度进行网格划分,因此计算量很大。

3.2 改进的震源扫描定位算法

王云宏[31]将震源扫描算法与DIRECT算法思想结合,提出基于DIRECT算法的微震快速网格震源搜索定位方法。该方法在求解震源扫描算法目标函数最优解时,不需要划分网格大小,直接搜索目标函数最优解。其基本步骤是:① 用式(11)构造目标函数;②根据实际资料,确定震源(x0,y0,z0)和发震时刻t0的大致区间范围;③ 利用DIRECT算法求解目标函数最优解,从而得到最终的定位结果。其特点是在远离最优解的地方搜索密度小,加快了计算效率;在靠近最优解的地方搜索密度大,提高了反演精度。该方法适合大空间内多事件、多参数的震源反演,不但为水力压裂微震监测实时、快速处理大量数据提供了可能,而且提高了定位精度。

何勇[32]等针对震源扫描算法对低信噪比监测资料定位效果不佳的问题,在震源扫描算法的“亮度”函数中引入了震源球的压缩和膨胀初动极性改正项,从而提出了一种改进的震源扫描定位方法。该方法将原来的只对震源位置(x0,y0,z0)和发震时刻t0进行扫描的算法,发展为一种对震源位置、发震时刻和震源机制同时进行扫描的算法,不但有效地提高了压制噪声和精确定位的能力,而且还给出了震源机制。

4 结语

地震定位作为微震监测技术研究的一项核心工作,长期以来诸多学者把加快数据处理、提高定位精度作为主要的研究方向。在线性定位法、非线性定位法等常规地震定位方法研究的基础上,发展了诸多适合微震定位的方法;其中无需预先测速的微震震源定位方法和震源扫描叠加算法,适合用于速度模型不准确或震相数据不清晰的条件下,比其他方法具有更高的定位精度。

随着采煤引发的冲击地压和煤层气开采过程中水力压裂效果监测的需要,对微震定位提出了更高要求。未来发展具有较强的抗噪性、消除速度模型等因素对地震定位的影响,将成为今后微震定位研究的主流方向之一。此外,为实现快速、准确的自动微震定位,还需要发展高质量的震相识别技术。

[1] 李楠.微震震源定位的关键因素作用机制及可靠性研究[D].徐州:中国矿业大学,2014.

[2] 刘高峰,翁红波,宋志敏,等.煤层气井水力压裂效果评价与消突时间预测研究[J].河南理工大学学报(自然科学版),2015,34(6):753-758.LIUGaofeng,WENGHongbo,SONGZhimin,et al.StudyonHydraulicFracturingEffectEvaluationandOutburstEliminationTimePredictioninCoalbedMethaneWells[J].JournalofHenanPolytechnicUniversity(NaturalScience),2015,34 (6):753-758.

[3] 董陇军,李夕兵,唐礼忠.影响微震震源定位精度的主要因素分析[J].科技导报,2013，31(24):26-32.DONGLongjun,LIXibing,TANGLizhong.MainInfluencingFactorsfortheAccuracyofMicroseismicSourceLocation[J].Science&TechnologyReview,2013,31(24):26-32.

[4] 李自红.临汾盆地地壳精细结构探测与孕震构造研究[D].太原:太原理工大学,2014.

[5] 孔韩东,边银菊,刘瑞丰,等.地震定位方法研究进展[J].地震地磁观测与研究,2017,38(4):81-92.KONGHandong,BIANYinju,LIURuifeng,et al.ReviewofSeismicLocationStudy[J].SeismologicalandGeomagneticObservationandResearch,2017,38(4):81-92.

[6]GEIGERL.ProbabilityMethodfortheDeterminationofEarthquakeEpicentersfromArrivalTimeOnly[J].BullSt.LouisUniv,1912,8(1):60-71.

[7]LEEWHK,LAHRJC.HYPO71 (revised):AComputerProgramforDeterminingHypocenter,Magnitude,andFirstMotionPatternofLocalEarthquakes[R].USGeolSurvOpen-FileRept,1975:1-116.

[8]EATONJP.HYPOLAYR:AComputerProgramforDeterminingHypocentersofLocalEarthquakesinanEarthConsistingofUniformFlatLayersoveraHalfSpace:R290-1253[R].US:DepartmentoftheInterior,GeologicalSurvey,NationalCenterforEarthquakeResearch,1969:1-117.

[9]LAHRJC.HYPOELLIPSE/version2.0:AComputerProgramforDeterminingLocalEarthquakeHypocentralParameters,Amgnitude,andFirstMotionPatterns[R].US:USGeologicalSurvey1989:1-119.

[10]KLEINFW.User’sGuidetoHYPOINVERSE-2000:AFortranProgramtoSolveforEarthquakeLocationsandMagnitudes[R].US:USGeologicalSurvey,2002.

[11]LIENERTBR,BERGE,FRAZERLN.HYPOCENTER:AnEarthquakeLocationMethodUsingCentered,Scaled,andAdaptivelyDampedLeastSquares[J].BulletinoftheSeismologicalSocietyofAmerica,1986,76(3):771-783.

[12]NELSONGD,VIDALEJE.EarthquakeLocationsby3-DFinite-DifferenceTravelTimes[J].BulletinoftheSeismologicalSocietyofAmerica,1990,80(2):395-410.

[13] 张少泉.近震分析[M].北京:地震出版社,1977.

[14] 杨文东,金星,李山有,等.地震定位研究及应用综述[J].地震工程与工程振动,2005,25(1):14-20.YANGWendong,JINXing,LIShanyou,et al.StudyofSeismicLocationMethods[J].EarthquakeEngineeringandEngineeringVibration,2005,25(1):14-20.

[15]DOUGLASA.JointEpicentreDetermination[J].Nature,1967,215(5096):47-48.

[16]DEWEYJW.SeismicityandTectonicsofWesternVenezuela[J].BulletinoftheSeismologicalSocietyofAmerica.1972,62(6):1711-1751.

[17]PAVLISGL,BOOKERJR.ProgressiveMultipleEventLocation(PMEL)[J].BulletinoftheSeismologicalSocietyofAmerica.1983,73(6A):1753-1777.

[18]PUJOLJ.JointEventLocation-TheJHDTechniqueandApplicationstoDatafromLocalSeismicNetworks[M].SpringerNetherlands:AdvancesinSeismicEventLocation,2000:163-204.

[19] 王椿镛,王溪莉,颜其中.昆明地震台网多事件定位问题的初步研究[J].地震学报,1993,15(2):136-145.

[20]CROSSONRS.CrustalStructureModelingofEarthquakeData:1.SimultaneousLeastSquaresEstimationofHypocenterandVelocityParameters[J].JournalofGeophysicalResearch.1976,81(17):3036-3046.

[21]AKIK,LEEW.DeterminationofThree-DimensionalVelocityAnomaliesunderaSeismicArrayUsingFirstPArrivalTimesfromLocalEarthquakes:1.AHomogeneousInitialModel[J].JournalofGeophysicalResearch.1976,81(23):4381-4399.

[22]KISSLINGE,ELLSWORTHWL,EBERHARTPD,et al.InitialReferenceModelsinLocalEarthquakeTomography[J].JournalofGeophysicalResearchSolidEarth,1994,99(B10):19635-19646.

[23]SPENCEW.RelativeEpicenterDeterminationUsingP-WaveArrival-TimeDifferences[J].BulletinoftheSeismologicalSocietyofAmerica,1980,70(1):171-183.

[24] 周仕勇,许忠淮,韩京,等.主地震定位法分析以及1997年新疆伽师强震群高精度定位[J].地震学报,1999,21(3):258-265.ZHOUShiyong,XUZhonghuai,HANJing,et al.AnalysisontheMasterEventMethodandPreciseLocationofthe1997JiashiStrongEarthquakeSwarmofWesternChina[J].ActaSeismologicaSinica,1999,21(3):258-265.

[25]WALDHAUSERF,ELLSWORTHWL.ADouble-DifferenceEarthquakeLocationAlgorithm:MethodandApplicationtotheNorthernHaywardFault,California[J].BulletinoftheSeismologicalSocietyofAmerica.2000,90(6):1353-1368.

[26]ZHANGH,THURBERCH.Double-differenceTomography:TheMethodandItsApplicationtotheHaywardFault,California[J].BulletinoftheSeismologicalSocietyofAmerica.2003,93(5):1875-1889.

[27]LIXB,DONGLJ.ComparisonofTwoMethodsinAcousticEmissionSourceLocationUsingFourSensorswithoutMeasuringSonicSpeed[J].SensorLetters.2011,9(5):2025-2029.

[28] 董陇军,李夕兵,唐礼忠,等. 无需预先测速的微震震源定位的数学形式及震源参数确定[J]. 岩石力学与工程学报,2011(10):2057-2067.DONGLongjun,LIXibing,TANGLizhong,et al.MathematicalFunctionsandPapametersforMicroseismicSourceLocationWithoutPre-measuringSpeed[J].ChineseJournalofRockMechanicsandEngineering,2011(10):2057-2067.

[29]KAOH,SHANSJ.TheSource-ScanningAlgorithm:MappingtheDistributionofSeismicSourcesinTimeandSpace[J].GeophysicalJournalInternational,2004,157(2):589-594.

[30] 李文军,陈棋福.用震源扫描算法(SSA)进行微震的定位[J].地震,2006,26(3):107-115.LIWenjun,CHENQifu.Micro-earthquakeLocationbyMeansofSource-scanningAlgorithm[J].Earthquake,2006,26(3):107-115

[31] 王云宏.基于DIRECT算法的微震震源快速网格搜索定位方法研究[J].地球物理学进展.2016,31(4):1700-1708.WANGYunhong.Grid-searchMethodonMicro-seismicSourceFastLocationBasedonDIRECTAlgorithm[J].ProgressinGeophysics,2016,31(4),1700-1708.

[32] 何勇,张建中,李同宇.改进的震源扫描定位方法[J].中国科技论文,2016(21):2450-2455.HEYong,ZHANGJianzhong,LITongyu.ImprovedSource-scanningAlgorithmforMicroseismicLocation[J].ChinaSciencepaper,2016,11(21):2450-2455.