摘要：利用Pn-Pg震相组合的局部搜索法，采用5种地壳速度模型对2013年吉林前郭5.8级地震震群进行震源深度测定。计算结果显示，5种结果的震源深度的主要分布特征相似。对选取的32个地震震源深度进行重新测定，结果最小为6 km，最大14 km，主要集中在10 km左右，系统深度偏差约1.74 km，震源深度分辨率约1.58km，重定位后地震震中优势分布呈北西向。与当地的沉积层厚度分布对比，这些地震最有可能发生在盆地基底的基岩顶层，而不在油气圈闭的沉积盖层内，并且地震震群分布走向与通榆一长春断裂走向相一致。结合该地区的历史地震活动与构造关系，最后推测这些地震与北西向的通榆一长春断裂构造运动密切相关，而不是油田作业的直接产物。关键词：震源深度；Pn - Pg算法；定位精度；油田作业；构造运动；前郭震群
中图分类号：P315.3 文献标识码：A 文章编号：1000-0666（2015）02-0211-10
0 引言
据《大金国志》记载，公元1119年2月19日～3月19日，在吉林省前郭卡拉木地区曾发生过一次强烈地震，经计算分析和综合考证，其震级为6%，烈度长轴为北西向（唐雅芝，1990）。另据周本刚和高名修（1992）不完全统计，在松辽盆地距离前郭县250 km范围内，历史上至少发生过5次M≥6的地震，分别是1119年前郭M6 3/4地震、1940年通辽M6.0地震、1941年绥化M.0地震、1942年通辽M6.0地震和1942年绥化M6.0地震。2006年3月31日吉林乾安一前郭还发生了一次M5.0地震，是1967～2006年吉林省发生的最大地震，通过75个现场宏观烈度考察点的结果，显示该次地震的烈度长轴方向也是北西向（盘晓东等，2007）。2013年10月31日吉林前郭发生了M5.5和M5.0地震后，于11月22、23日又发生了M.3、M5.8和M5.0地震，期间还有不少小震事件。地震发生后，一些科研院所在应急T作阶段给出了部分巾强震的震源深度，初步结果为5 km左右，通过这个深度特征认为这些地震可能与该地区的采油作业有关。由于前郭地区位于松辽盆地中部，在第二松花江断裂（图1）北侧就是大庆油田采油区。而一些地震勘探资料显示，乾安地区的采油井深度多在3800 m左右（高金哲等，2013）。这些地震是否归因于采油作业的影响，其深度定位的绝对值与精度对结论有决定性作用。由于松辽盆地第四纪沉积比较薄，厚约75～200 m（裘善文等，2012），而作为我国中新生代陆相人型沉积盆地，沉积面积居第三位，其中白垩系沉积厚度人于7000 m（杨令宾，1995），部分测线显示最大沉积平均厚度6093 m，考虑地表剥蚀效应后求得沉积物最大平均总厚度及基底总沉积量为7928. 58 m（胡望水等，2005），因此过浅的地震震源深度是值得怀疑的。
本文将利用地壳内的直达波震相Pg和沿Moho面滑移的首波Pn进行震相组合，重新测定2013年前郭震群进行初始破裂点深度。这种Pn-Pg组合方法发挥了这两个震相的离源射线上下反向而对震源深度变化更为敏感，以及初至震相到时提取精度高的双重优点。本文进一步利用初始定位的震巾位置，在震中附近采用局部搜索法，以获取Pn-Pg最小走时残差的震源深度和震中位置，并对不同速度模型的结果进行讨论，以此给出前郭震群比较合理的震源深度分布特征。
1 区域构造特征及Pn-Pg局部搜索法
1.1 前郭地区局部地质构造
前郭地区位于松辽盆地中部，松辽盆地产生于燕山运动时期，从晚侏罗世到新生代，盆地经历了裂谷型断陷、大型拗陷和萎缩上升阶段，它是一个人型中新生代陆相沉积盆地（胡望水等，2005）。喜山运动的三次运动幕在松辽盆地都有表现。第一、三幕运动主要受印度板块影响，产生的区域主压应力为NE-SW向，第二幕运动主要受太平洋板块影响，产生的区域主压应力方向为NW-SE或近东西向。在卫星遥感图像上可以地看到，盆地在喜山第二运动幕中产生了系列轻微的NW向构造形迹线。而全新世以来，盆地内的江河流向与湖沼分布都受到了NW向构造线性体的控制，如人庆油田及邻区地表的湖泡展布，这些线性体时隐时现地切过了大庆长垣表层第四系。在扶余幅卫星图像上，不仅可以看到乾安周缘大火小小星罗棋布的湖泡呈NW向串珠状排列，还可以看到靠近农安附近呈NW向展布的波罗泡群及其一带的全新统Q4沉积，切割并覆盖了盆地东南隆起区中的NE向构造群。在通榆幅卫片上，还可以解译出NW向的构造线性体，斜穿了NE向弧形构造带，可见其形成时间最晚（祁林等，1992）。
对现今或喜山三幕区域应力场表现最为敏感和直观的就是近代地震。垂直于NE-SW区域主压应力作用方向发育的NW向构造线和断裂面，是地应力不断积累的地带。刘天革和隋红（ 1985）统计东北地区历史地震活动的研究结果得出，整个东北地区共记载有60多次M≥4地震，其中32次是沿着NW方向迁移发震或成对出现。其余20多次则发生在NW向断裂与NE向断裂的交切点上。这些强震的特点是余震分布往往呈NW向展布，等震线形状往往呈NE、NW双轴型。显然，东北地区发震机制是山于NW向断裂应力积累而活动，或者是由此诱发了NE向断裂共同参与活动而发震。而位于松辽盆地南部地区的乾安油田，油田周围所有的泡子都逐渐下洞或消亡，说明该NW向新构造圈闭仍处于隆起状态中（祁林等，1992）。
在本次地震震群地区，分布着通榆长春断裂、双山一前郭断裂、人安一都德断裂、第二松花江断裂以及其他邻近断裂体系（图1）。图1中虚线框里的断裂分布抽取自盘晓东等（2007）。其中通榆长春断裂为一条隐伏断裂，断裂切过Mo-ho面，是控制中央凹陷带南部长岭凹陷的两个次级凹陷（即乾安次凹陷和黑帝庙次凹陷）的分界断裂。由于第四纪以来的多次活动，导致沿通榆一长岭怀德一线出现了一个活动构造地貌单元，即松辽分水岭，形成了北部区的松花江水系与南部区的辽河水系。而大安都德断裂沿NNE向从吉林省的人安一带向北到黑龙江省的都德一带，主要活动时期为中侏罗世到白垩纪及晚更新世末期，全新世以来也表现出较强的活动性。2005年7月2s日发生于黑龙江省林甸5.1级地震就与该断裂的活动有关。2006年3月31日乾安一前郭5.0级地震则发生在这两条断裂的交汇部位（盘晓东等，2007）。 1.2 Pn-Pg局部搜索法
变化比值，此处称之为灵敏度k，其值与震中距△，震源深度h有关。
本文利用该方法进行深度定位。考虑吉林省地震局人员在应急编目过程中，因时间仓促等因素可能影响震相拾取精度，以及地震震中定位的误差影响，本文对初至信噪比较高的附录表中的地震逐一进行了Pg、Pn震相重新识别，并结合100 km以内的直达剪切波Sg震相，利用中国地震局推广使用的MSDP软件中的单纯形法进行重新定位（速度模型为华南走时表），并将该定位震巾结果作为搜索参考点。搜索过程取参考点水平面±0.2。范围，步长0.01°，深度取0～30 km，步长1 km，即52111个三维网格点。计算每个网格点到每个台站的Pg、Pn理论走时，并与台站人工提取的Pg、Pn走时进行差值比较，单个网格点的走时残差定义如下：
将每个网格点的走时进行对比，可以获取到最小走时残差的网格点的走时、经纬度和深度值。虽然经纬度精度只有2位小数点，深度值只精确到整数，但已是本方法的精度限，而这个精度已能为本文的后续分析给予足够好的支持。对于2013年12月8H前郭Ms3.1地震，其分析所用到的流动台站和遥测台网台站分布如图1所示。沿图1中震巾位置两侧及垂向方向进行局部网格化，计算出各节点的走时残差，其切片如图2a所示。可见其走时残差分布呈近似椭球状，椭球中心点残差最小。但受速度模型分层（本例采用模型d）的影响，椭球是分层连续的。
2 前郭震群的地震深度测定及其精度分析
2.1 前郭震群的地震深度测定
本文处理了2013年10月31～12月27日前郭地震序列波形资料，由于方法的限制，只拾取至少有一个Pn震相清楚的地震地震事件，总共获取了32个满足条件的地震事件（表2）。极少量台站（如SYZT台）存在GPS时钟对时不准的问题。对32个满足条件的地震沿震中距作Pg、Pn震相时距散点图（图3），显示绝人多数没有显著的GPS时钟差异问题，所以只剔除了个别有问题台站的震相，走时点状图不完全符合和达理论走时曲线，其原因可能与地壳各向异性有关。
在地震定位结果中，地壳速度模型会对结果有较大的影响。由于大庆地区石油物探剖面较多，可以为本文提供参考。刘洋等（2008）给出的大庆地区上地壳P波速度显示，在地表3 km以下，速度剖面大多数超过4.5 km/s。松辽盆地及其南北延展地区是一个Moho面拱起带，它由3个上地幔隆起组成，走向北东20°左右，地壳厚29～33 km，其东西两面的地壳厚度在35 km以上（程学儒，1982）。松辽盆地东北隆起区Moho界面双程反射旅行时变化范围为9.6～11.0 s，深度范围为30～34 km（杨宝俊等，2003）。在松辽平原形成过程中，白垩系沉积物起着重要作用。如果以各地质时代地层厚度占中新生代地层总厚度的百分比作为沉积物在平原形成过程中的贡献率，则白垩系沉积物贡献率达86.37%，第三系为7.85%，第四系为5.78%（杨令宾，1995），这种成分上的比例也为速度模型提供约束。东部张广才岭段的地壳厚于盆地段，为34～37 km，上地壳厚14～15 km，以中酸性与酸性岩石为主，Vp为6.02 km/s；中地壳较薄5.5～7 km，以中酸性成分的变质岩为主，Vp为6.5 km/s；下地壳11N15 km，以基性成分的变质岩为主，Vp为6.85 km/s。盆地内盖层结构较简单，速度随深度递增，4.75～5.25 km/s，Moho面速度8.00 km/s（傅维洲等，1998）。对比华南广东地区与松辽盆地的Moho面深度分布（滕古文，曾融生，2002），发现二者相差较小，松辽盆地略深约1 km。因此，在参考华南走时表模型的基础上，给出了另外5种地壳速度分层模型（表2、图4），其中模型b、d、e含有地表低速层。
通过6种速度模型（表3巾计算32个地震震源深度），模型a、模型d、模型f的结果见表3。其他3种情况的结果，考虑篇幅的限制，仪将结果绘于图5、6。图5、6显示：①MSDP的单纯形法给出的结果相对偏浅；②模型加入地表低速层后，深度值系统性偏浅约1 km；③Moho面加深后，深度值系统性偏深约1 km；④增加了地表低速层后，人多数情况可降低走时残差，特别是7个流动台震相加入定位的情况更为明显（表3中加圆括号的ID号事件）。但如果进一步增加盖层厚度并降低盖层的P波速度，则残差变大。
不同速度模型的计算结果反映了一些共性特征：在前一组5级地震后和后一组5级地震后，小震深度都有变浅的趋势。对于ID为10号的Ms2.3地震，各种模型计算的结果都是最小的，范围6～8 km，最小值为6 km。这个深度为中生代沉积的底部、基岩顶层深度。因此利用Pn-Pg组合定位的局部搜索方法，给出的32个地震深度都在6 km以下，位于盆地巾生代沉积的底部、基岩顶层交界及以下。另外从搜索的震中分布显示，5种模型给出的32个地震震中具有北西向分布特征。人震级的地震偏南深度偏深，小震级的地震偏北偏浅。
2.2 前郭震群的地震深度精度分析
笔者对流动台的Pg到时数据进行先后排序，发现流动台L2202离多个地震的震中非常近，约3 km以内。为此，挑出了震中位置几乎相同的3组地震，ID号分别为10和16、25和26、27和29进行单台Pg-Sg走时对比。其中第1组前者深度较小，第2组二者深度相同，第3组后者深度较小。由于采样率为100 sps，每个地震到L2202流动台的直达波Pg和Sg震相经过局部放大可精确提取初至时刻点，时间拾取精度可以达到0.01 s。在实际拾取过程中，对比前后波形差异，取波形起始较尖锐并且信噪比人于3的点作为震相的初至点。图7中黑实线为各组地震中的前1个地震波形，红虚线为后1个地震波形，波形统一将Pg初至调整到第100点进行对齐，每组地震的Sg-Pg到时差见图7a'～c' ，相应的波形见图7a～c，其中括号内数字为深度，前面数字为Sg-Pg的到时差。 对比图7a -a'可以发现，7 km深度的地震（ID=10）震中距较远，但Sg-Pg到时差小，波形巾Sg提前。对比图7b-b'，两个地震震源深度相当，震中距远1 km的（ID=25）地震Sg-Pg到时差略大，两次地震记录的波形Sg几乎同时到达。对比图7c-c'，震巾距小的地震（ID=27）震源深度深，Sg-Pg到时差略大，波形中Sg延后。由于同震巾点的地震，震源深度没有明显差异，此处未能得到比较。但是从上述三组事件可以较清晰地看出，用Pn-Pg搜索的深度与极近台Sg-Pg到时差的对比结论是一致的。因此本文给出的各模型结果即使相互间可能有一定的系统偏差，但不同模型确定的地震的相对深度是比较可信的。3组地震中，最大深度12 km，最小7 km，进一步考虑震中距不同的影响，这5 km的深度差异识别度应仅是下限。对于第三组的3 km的深度差，地震波分析也给出了分辨效果。因此综合震巾距的影响，其分辨率约为0.5×[1/3（12-7）+1/3（11-8）]=1.58 km。而对于系统偏差，参照图5b可以发现其残差均值在0.4 s附近，由于这是Pn、Pg反向射线走时的差异，所以单Pg残差上限可估计为（1-0.333）×0.4 s=0.28 s，若上地壳P波取6.2 km/s速度，其系统深度误差为1.74 km。而利用1.2节的速度模型，可以推出其虚波速度为8.232 km/s，对于ID号为16、26的地震，按虚波速度计算的深度结果与本文是基本一致的（剔除第四纪盖层厚度后吻合度更好）。
3 认识与结论
虽有研究显示，油气开采及其注水增产作业可能会诱发小地震（刘元生，王六桥，1985），但Davis等（1995）提出判别流体抽取导致地震的发生需考虑几条准则，其巾3条是：
（1）是否有些地震发生在抽取流体的深度上；
（2）是否是该地区第一次大地震；
（3）是否在抽取流体后发生。
因此在石油开采过程具体资料比较欠缺的情况下，震源深度测定及测定精度对判断地震的主要发震因素非常重要。
本文利用空间局部网格点有穷寻优的Pn-Pg组合法对2013年10月至12月吉林前郭震群进行了深度测定，并利用流动台的观测波形对测定结果进行检验，给出了测定的系统误差与单次地震的深度误差确定性估计。计算结果表明，这32个地震的深度在6 km及以下，超过松辽盆地巾生代沉积的底部、基岩顶层深度，特别是Ms≥3.3地震的深度超过8 km。因此这些地震最有可能发生在盆地基底的基岩顶层，而不在油气圈闭的沉积盖层内。由于地震震群分布人致呈北西向，与通榆一长春断裂这条切过Moho面的隐伏断裂走向相一致，因此本文推测这些地震与北西向构造活动有关。
另外，据中国地震局地球物理研究所（2013a，b）利用波形拟合方法给出的震源机制解，10月31日5.5级地震的震源机制解为：节面I（strike330°，dip60°，slip50°），节面 II（ strike209°，dip48°，slip138°）；11月23日5.8级地震的震源机制解为：节面I （strike200°，dip70°，slip140°），节面II（strike3060， dip53°， slip25°）。如果按本文北西向的震巾分布作为参考，则上述第一个地震的主破裂面为节面I，第二个地震的主破裂面为节面II，都带有一定的逆冲分量，西南盘向东向走滑并兼带仰冲。这与该区域全新世NW向新构造圈闭仍处于隆起状态的结论吻合。
由于前郭地区1119年曾发生6 3/4级地震，2003年4月16日至7月31日吉林松原出现震群活动后，2005年7月25日在大安一都德断裂的林甸发生5.1级地震，其后2006年3月31日乾安一前郭发生5.0级地震，以及2013年10至12月的多个5级地震的发生，似乎都是北北东向大安一都德断裂和北西向通榆一长春断裂构造活动的结果（盘晓东等，2007）。因此综合前郭地区的历史地震、本次震群的震源深度、震中分布及局部构造运动特征，笔者认为发生在2013年10至12月的前郭震群是构造运动的产物，油田作业不是其主要发震因素。