• / 50
  • 下载费用:1 下载币  

第八章 地震定位与地球速度反演3

关 键 词:
物探 地震资料解释 地震处理 反演
资源描述:
89第八章 走时数据的反演 ...................................................................................................................... 球层介质的速度反演 ................................................................................................................点法(古登堡法) ..................................................................................................... 法 ............................................................................................... 单台地震定位及 定位 ....................................................................................台定位法 .......................................................................................................................站求取震中位置 ........................................................................................................用三台站求解地震位置的计算机算法 ........................................................... 地震定位结果与台站分布的关系 ....................................................................................... 迭代定位方法 ..............................................................................................................................代定位算法的基本思路 .........................................................................................代方法的具体实现 .................................................................................................... 相对定位方法 ............................................................................................................................ 章 走 时 数 据 的 反 演在前两章我们根据已知的速度结构对射线追踪和走时曲线的计算问题进行探讨,推导了波速只随深度变化的一维速度模型的射线追踪的表达式。三维结构的射线追踪虽然较复杂,但它遵循的原则是类似的。现在我们来研究根据观测到时数据得到速度结构和地震位置的问题。前面两章是已知速度结构和地震位置,如何计算地震波走时的问题,而本章是根据地震波走时得到地震结构和震源位置的问题,是前面两章的反问题。在这个问题的研究中,地震学家通常把问题进行简化:第一种简化是对于震源位置已经清楚的地震,根据地震波走时数据求解速度结构,第二种简化是已经得出了该地区的速度结构,根据地震波到时求解地震位置。限于篇幅,对于第一种简化,本书仅讲授球层介质一维速度结构的反演。对于第二种简化,我们讲授地震定位的基本方法。 点 法 ( 古 登 堡 法 )这个方法是利用走时曲线的拐点处与从震源处水平方向射出的射线相对应,据此可求出震源处的波速。根据第七章球层介质中的 律有:\* 8h 代表震源处的参数, 都是常数,从震源向不同方向射,出的射线, 值不同,相应的射线参数 p 也不同,当 =90时, 达到极大值,hi 因此有* 82律(7, ,有 82应于走时曲线的拐点。因此走时曲线的拐点与从震源处水平射出的2样,我们只要求得某地震的震源深度及有观察分析归纳得到其走时曲线,找出走时曲线的拐点,并确定该点的走时曲线的斜率 ,则\* 8 点有 ,因此有2 8中,h 为震源深度, 为走时曲线拐点 M 点的视速度(参看第七章球层介质的律),R 为地球半径。此方法原理清楚、方法简单、计算方便。但由于地震发生的深度大约在 0用这种方法只能求出 0的波速。由于拐点不易找准,得到的精度较差,它又要求对应每个深度的地震就要总结出一条走时曲线,因而资料分析工作量很大。 法考虑球面介质,速度随半径 r 变化[v( r )],根据 中得出震中距的表达式有\* 128中 = , 为射线转折点距地心的距离,p= ,R 为地球半径。改变积分变数,( \* 02128中 。0在此介绍一积分式,, \* 01218为转折半径为 的射线参数,所以\* 8表示11 至到转折半径为 的所有射线的积分。将\* 81作用于震中距表达式\* 8的两边,得到\*  0101 2122 p 8程的左边变为 ,考虑到 ,可以用分01 122111''* (8* 8的第一项代入积分上下限,由于 )=0,且在 p=,半径为 R,Δ=0,因此第一项消失。仅剩下第二项,可写成:0\* 101)(8* 8右侧为一双重积分式,改变积分变数,得到:1111 111 1111221 222 2241 112()l|pr r     1\* 8里采用了不定积分公式:   01 1()8\* 8结合,得出\* 10118 ,可以得到:101\* 1就是射线最深点的半径的表达式。此表达式可以根据观测走时曲线 t( )得出。如果 t( )为一平滑曲线,我们可由其切线斜率 得到 p( )(由于 p= )。() 8中的 ,为在震中距 的走时曲线斜11p1 的积分式由震中距 Δ 从 0 到 Δ 1 的 p 值,得到 ,根据\* 8求出 ,即其最深转折点。根据转折点为 的 ,r 1得出其对应此深度的速度\* 8依此类推,得出速度与深度的分布图。下面归纳一下求速度分布的具体步骤1、 做出 0- 范围内的走时曲线,即 线,震源要取在地表上,否则要进行修正;在该步骤中,通常在计算机上利用三次自然样条函数对一元函数进行成组插值及微分,从而使得走时曲线具有较好的平滑性能,采用该方法能保证所插值的函数及其一阶导数、二阶导数连续。2、 由 线求曲线上各点的斜率,从而做出 曲线;该步骤求出震中距的分布。3、 在 0- 范围内任取 ,取出 点对应的 ;这就是111t。11)(、 将 曲线中 0- 段除以常数 即 ,得到即 曲线,从而求出 曲线;11p1震中距的变化。15、 在 0- 范围内求出 曲线下的面积,即求11p,这就是 ;在计算积分的10101过程中,通常采用 分数值公式进行。根据\* 8求出 ,然后根据\* r(8得到 所对应的速度 。1取不同的 值,重复上述步骤,从而求出不同的 值所对应的速度值 ;1用的数据是北京台网、大连、长春、牡丹江等 26 个台站的记录,并选用从北京地区到阿留申群岛西端的 300 多个地震,经过震源深度及剥壳校正后(地球平均厚度 35公里)。得到的北京 波走时曲线。数据如表 8北京- 萨哈林剖面P 波走时观测值。表 8京- 萨哈林剖面 P 波走时观测值∆(°) t(s) ∆(°) t(s) 法,利用 的三次样条插值和辛普森积分方法,编制成 序。程序如下: %调用数据x0=,1)';%震中距y0=,2)';%走时.2;x=0.2:7.8;n=x);[y,x0,y0,x); %对数值进行样条插值,得到插值点的 阶导数 阶导数 ) %x,y,'r*');实际观测曲线') 震中距 ');走时');R=6371;r1=,(:2:n))); %射线能够达到的最低点的半径v=,(:2:n))); %射线最低点的速度:2: %视为常数 f=: *(f(1)+f(2*f(3:2:+4*f(2:2:)*80; %采用 式进行积分  )(2(4)6)( R/nt/ %采用公式求解射线最深点的半径)1Srv( %采用球层 律得到最深点的速度; 函数的作用是:将角度转换为弧度 ,利用公式 1; %序号+1%v,.') %绘制速度和深度剖面图' %将 y 轴的方向进行翻转利用上述程序调用表 8的北京- 萨哈林地区的 P 波走时曲线的数据。对观测曲线进行上述程序来反演即可得到速度剖面如图 8京- 萨哈林地区 P 波走时观测反演速度剖面。8京- 萨哈林地区 P 波走时观测反演速度剖面对反演数据进行重新插值处理,选择出合适的数据,分别选择如下的数据(20,45,85,125,145,165,265,365,405,445,495,565,705,875,965,1025,1165在上述 数后加入如:20,45,85,125,145,165,265,365,405,445,495,565,705,875,965,1025,1165];[y,v,x);[v'][x',y']并利用得到的数据进行绘图得到图 8点插值的反演数据图图8些公式是有问题的。在这种情况下, 不连续,没有)(p唯一的解。对这种情况,传统的做法是设想在低速区存在若干有不同速度的均匀层。由于没有射线在这些层里转折,因此低速区的这些层可以任意移来移去,穿过低速区的射线的积分的走时和震中距不会改变。不论解析多么完全,在地震学中很少用赫格劳兹—维歇尔(式,这至少有两个原因:首先,定我们有已知的连续的 曲线,而实际上我()t们总是只有有限的走时点。这意味着走时曲线必须在这些数据之间进行内插,这些内插方案的不同导致于有不同的速度剖面。实际上,会有无数个稍有不同的速度模型,它们都适合于有限数目的 点。然而,更严重的问题是实际的()t地震数据一般多少都有噪声,自身相互矛盾。图 8左边的例子中,小的时间变动导致 点的波动,不可能把这些点()t与能得到期望的 1度模型的光滑的、物理上可靠的 曲线相联系。在右()们注意到可以把几次地震的数据结合在一起,用这些数据确定“平均”的速度剖面。图 8呈现离散的走时观测结果往往使速度剖面的反演复杂化由于这两个原因,式不能直接应用。下面我们对地震学家反演这些不完善的数据的一些方法进行讨论。式的主要用途是阐明精确的 曲()此,许多反演的策略是把寻找最佳速度模型的问题简化为寻找最佳的 曲线的问题。()t台地震定位及 定位根据走时数据确定地震的位置是地震学中的一个“古老”的,富有挑战性的问题,一直是地震学研究的一个重要方面。地震通常按发震时间和震源位置来定义。震源是地震的位置 ,而震中是震源正上方地表上的一个点),(地震定位中,通常把地震作为点源来对待。对于破裂几十到几百公里),(源不一定是地震的“中心”,宁可把它看成是地震发生时,地震能量开始辐射的那个点。由于破裂的速度小于 波速度,可以不管地震最终的根据初至到时来确定震源。标准目录给出的地震资料,例如震中的初步确定(依据的是高频体波的走时,不能把这些发震时间和震源位置与长周期反演结果相混淆,后者往往给出地震的矩心时间和位置,节我们首先讨论单台定位法和多台定位的 法。 台 定 位 法单台定位的原理就是根据纵波初动确定出震中方位角,然后根据震相到时(走时表等)确定出震源到台站的距离,根据水平向和垂直向记录的数值,可以求得视入射角,从而根据地震波速度估计出真入射角。根据震源到台站的距离和真入射角求得地震的震中距和震源深度。注意震中方位角是指通过台站子午线与地震台站到震中连线间的夹角,沿顺时针量取为正。由图 8示,设地震台站的 P 波的北向初动半周期位移为 ,东向初动半则水平向位移矢量与北向的夹角可以由下式来计算 8在 可以用 uE,求得,因为 数根据的正负号给出到 的范围,将其转换为度,并且对于小于 0 度的角度加上 360 度,就转换为0~360 度范围。知道了水平向矢量方向,地震震中应该位于该矢量方向上,但究竟是沿着矢量方向还是与矢量方向相反?这需要采用垂直向记录。我们定义垂直向记录 向上为正。 垂直方向的初动半周期决定地震射线的方向。震射线的方向背向震中;当垂直初动向下时,地震波射线的方向指向震中。图 8一断层错动图,台 1 观测到的是压缩波(+),其 P 波垂直向记录向上,即射线是“离源”(即背向震中)运动;此时震中位于水平矢量的相反方向上。台 2 观测到的是膨胀波(-),其 P 波垂直向初动向下,即波射线为“向源”(指向震中)的运动。此时震中位于水平矢量方向上。根据这种分析,我们可以得到结论:垂直方向初动向上的台站,震中位于水平矢量方向的相反方向上;反之垂直方向初动向上的台站,震中位于水平矢量方向上。个台站接收到地震 P 波初动符号,判断地震所在方位的示意图第二步,确定震源到台站的距离设 别为该台站 S 波及 P 波的到时,则地震到达该台站的 S 波和 P 波的走时之差等于 S 波及 P 波的到时差,则有:\* 8通常定义\* 18物理量具有速度的量纲,可以根据当地的平均地壳 S 波速度和 P 波速度估计该参数。通常被称为虚波速度。设若 5.0 km/s,= 3.0 km/s 则虚波速度 k=s。根据台站上读取的该地震的 S 波和 P 波到时差和该地区的虚波速度,可求得震源到台站的距离 r,公式为…………………………………………………..\* 8震记录的垂直向和水平向的记录,求得视入射角,从而根据地球均匀速度模型估计真入射角。地震在该台站的水平向位移由南北向和东西向位移记录合成:\* 28视入射角为:\* 8据第四章的真入射角和视入射角的转换公式(4即 得到真入射角。第四步,根据真入射角和震源到台站的距离 r 求得震中距 和震源深度 h。\* 8* 8中方位角及震中距求得以后,便可决定出震中位置。这样应用一个台站的三分向记录就可以估计地震的大概位置了。在计算机高速发展的今天,计算机可以代替人工操作。我们通常给出台站的经纬度,需要求解地震的经纬度和深度。在近震的情况下,通常采用直角坐标来解决问题,因此需要进行经纬度转换为直角坐标系的转换。将观测点的经纬度 转换为直角坐标 有很多成熟的公式可以使用i,如,朱介寿等,1988)。由于本研究的区域范围较小,我们采用下述简单公式( 990)。将观测点的经纬度 换算为直角坐标 的公式如下:i,* 2008中 为直角坐标的原点。 , 以度为单位, xi,0,0i序如下:x,y]=根据坐标系原点地理坐标(换为直角坐标x,y%输入 输出向),向) %1度弧长y= %北向坐标(8第一式2; %纬度的平均值作为计算经度的大圆的半径x= %经向坐标(8第二式需将直角坐标系中的坐标 转换为经纬度 ,公式为:i,\* 00028中的符号及意义与上式相同。采用 行转换的程序如下:x,y,根据坐标系原点直角坐标x,入 x, 输出 %y/c1+ %(x,y)点的纬度,(8第一式2; %纬度的平均值作为计算经度的大圆的半径x/c1+ %(x,y)点的经度,(8第二式们以北纬 25°,东经 120°作为坐标原点,求北纬 26°,东经121°的直角坐标,运行“[x,y]=6,121,25,120)” ,则可以得到 x= y=结果,然后用此结果求解其经纬度,运行“[x,y,25,120)”的语句,可以得到北纬 26°,东经 121°的结果。读者也可以采用其他的数据进行检验。有了上面的坐标转换程序,我们可以编写根据单台记录求解地震位置的计算机程序为:5; 15; %台站的经纬度;;2; %地震台记录的南北分向、东西分向和上下分向;; %位为秒;; %平均的位为km/vp*( %虚波速度ue, %按(8求解方位角if(00) %如果地震深度超出可信范围,则强制给定地震深度) = 12;) = %估计地震发震时刻的平均值50 0 50 100 100 100 50 0 0 ]; %台站位置的50 100 100 100 50 0 0 0 50]; %台站位置的,9); %台站位置的此设置为零 %各个台站的直达 %各个台站的直达,3,g, %调用和达法求解地震位置在上面的程序中,我们均采用直角坐标。如果已知的台站坐标的经纬度,求解地震的经纬度,则需要采用(8首先转换为直角坐标,待求得地震位置后再采用(8转换为地震的经纬度地理坐标。震定位结果与台站分布的关系地震定位问题需要求解地震的空间位置和发震时刻。我们把这些参数称为模型,定义模型矢量:\* (),,(4321 8震相观测到时。为了得到参数 ,我们首先一维平层介质模型。这样,按照第 6、7 章的知识,根据 值和台站位置可以计算第 个台站的预测到时: 8里 是算子,其复杂程度由假定地球模型的复杂程度和地震波传播的射线路测到时与预测的理论到时之差为:\* 8里 是第 个台站的残差。在某种意义上来说,我们希望得到使观测到时与预注意, 是地球模型和每个台站位置的函数,根据,7 章的知识就可以精确计算。然而,走时与地震位置参数的非线关系使反演最佳地震模型的工作大大复杂化。这种非线性的影响基至在均匀速度进行定位的简单情况下,也是明显的,在这种情况下,座标为 的台站记录震源(,)震时刻为 T 的到时为:(,) 2()()(8里 是速度。在此方程中,到时 不论是与 还是 都不成线性关系。该结果们不能用解线性方程的方法来得到震源位置。现在有了功能强大的计算机,我们可在所有可能的地震位置和发震时间范围里作网格搜索,计算每个台站预测的到时。那么,我们可以求得使预测的时间 与观测的时间 最一致的特定的 。怎样规定“最”一致呢?流行的选择求使下式达最小的 m 值: 8里 是台站数目。把平均的平方残差 叫做方差。在数理统计上更普遍的式是用 来定义方差。这里 是自由度的数目( 为 减去拟合中的2df对通常需要解决的问题,拟合参数的数目比数据的数目少得多,所以 和 近似相等。在数理统计中 分布于自由度个数及置信度的关系如表 8示。表 8布与自由度和置信度的关系25%)2(50%)2(5%)25 是因为在求最小值的问题中,它使方程有简单的解析形式。如果 与 之间的拟合差是由 中非相关的随机噪声所么最小二乘法将有助于给出正确的答案。在本门课程的讨论中,如不特殊说明,我们假定误差是由服从高斯正态分布的随机噪声导致的。现在我们以一组台站,假设为均匀的地壳模型,以 1 网格点搜索地震的位置。序如下:0;9; %假定的震源位置,发震时刻假定为零S1%0;0; %假定的震源位置,发震时刻假定为零S2%0;0; %假定的震源位置,发震时刻假定为零; %假定的地震波速度%假定的台站位置N= %台站个数Pt=,1)^2+(,2)^2)./1,1))* %计算各台站的预测走时, %该问题需要求解地震的发震时刻和x,此有三个未知数 %给出最小值是一个较大的数;;t); %初始的位置x=0:1:100y=0:1:100x*,1),1)).^2+(y*,1),2)).^2)/ %按(8求得预测到时.1:^2)/ %按照(8计算平均残差if( ) = 进行:循环迭代,到导值和各台残差,方位角和离源角) [p,g,%计算各台定位残差 = %保存循环中每次计算的残差,= 保存最小残差值及对应的震源参数 %每次循环都要判断1次,显得太繁琐) = %)/()/ 1);) = )/()/) = )/()/ % 0;) = )/()/ %0;F%如经纬度修正值,大于系统设定值 ) > )% );% );% ) > ) = ) = )/()/) = )/()/) = )/()/) > ) = )/()/) = ) = )/()/) = )/()/irLo
展开阅读全文
  石油文库所有资源均是用户自行上传分享,仅供网友学习交流,未经上传用户书面授权,请勿作他用。
0条评论

还可以输入200字符

暂无评论,赶快抢占沙发吧。

关于本文
本文标题:第八章 地震定位与地球速度反演3
链接地址:http://oilwenku.com/p-63151.html
关于我们 - 网站声明 - 网站地图 - 资源地图 - 友情链接 - 网站客服客服 - 联系我们
copyright@ 2016-2020 石油文库网站版权所有
经营许可证编号:川B2-20120048,ICP备案号:蜀ICP备11026253号-10号
收起
展开