专利名称：一种用于地震勘探中射线追踪的方法
技术领域：
本发明属于石油地震勘探资料处理方法技术领域，尤其涉及一种最小旅行时地震射线追踪的方法。
背景技术：
地震勘探法，即按一定规律和方式排布的各地震波检测点接收由按一定规律或方式排布和移动的各放炮点在地表激发所产生的地震波的折射波和/或反射波，通过对检测到的由各地质构造层界面处对该地震波产生的折射波和/或反射波返回地表的时间和波形进行处理和分析，以便能准确了解和确定地质构造的形态和位置，是目前在石油和天然气等勘探中指导勘探钻井孔位部署的一种常用和有效的方法。由于地下各地质构造层在组成及密度、分布的均匀程度等诸方面的不同，放炮点激发产生的地震波在其中的传播速度相应也各不相同，各检测点检测到的由各地质构造层界面处产生而返回地表的相应折射波和/或反射波的时间早晚及波形等也相应有所差别。
地震波在地下传播的过程中，波穿过的路径和介质的速度决定了波的旅行时间。根据波动理论，波从激发点(震源点，如地震时的震中)出发，以介质速度在地下传播，经过一段时间，到达接收点(记录站点)，所走的路径总是旅行时间最短。这种最短路径(射线旅行时最短的路径)就是波的射线路径。如何从激发点到接收点确定波传播的路径和旅行时，就是射线追踪研究的主要内容。在地震偏移方法、层析成像和反演方法中，广泛采用射线追踪方法确定地震波旅行时，因此射线追踪方法是一个重要研究领域。现有技术中存在几种射线追踪的方法，如两点射线追踪法(Feregra，1980)，有限差分方法解程函方程法(Vidale，1988)，波前重建估算旅行时和振幅(Vinje等，1993)，波前追踪法(Bube等，2001)。
但，所有射线追踪技术主要面临的两个问题是1、射线在追踪过程中产生的盲区(即射线路径走不到的区域)和并非时间最短路径。射线的盲区主要发生在地层的低速度带，当射线遇到低速层时，绕过低速层，从高速度层穿过，因而产生了盲区。
2、有些射线追踪方法有可能找不到最短路径(旅行时间最小)，如四方网格打靶法。运算速度也是实际应用中必须考虑的问题。耗时过多的算法，尽管能克服射线的盲区问题，也无法在实际应用中得到广泛应用。
为了克服上述存在的问题，现有技术中提出了一种矩形网格三点扰动(中间点做水平位置的滑动)射线追踪法，该方法是针对速度模型离散的情况提出的，不论什么样的速度场分布或低速带的存在，该方法总能找到最短时间的路径，避免了射线盲区和追踪路径并非时间最短路径问题。其具体操作方法是1、对追踪区域作矩形网格化，每个网格上分布速度值。在震源点到接收点用直线联结，得到初始射线路径，初始路径与网格的交点分别为p1、p2、…pN。计算路径p1p2…pN总旅行时。
2、取初始路径与矩形网格三个交点pi-1，pi和pi+1，固定两个端点pi-1和pi+1，对中间点pi修正，给出一个扰动量(扰动的增量)对其扰动，计算相应的旅行时间。所有可能是最小旅行时的位置都进行这样的位置滑动，计算位置滑动后的旅行时，然后比较旅行时，确定最小旅行时，从而确定具有最小旅行时的位置qi，用qi代替pi的位置。
3、取qi，pi+1和pi+2三个点，固定两个端点qi和pi+2，重复2的扰动和旅行时计算。确定qi+1代替pi+1的位置。这样依次取三点修正直到修正到pN为止，并计算新路径p1q1q2…pN的旅行时T。
4、把p1q1q2…pN的连线当作新的初始射线路径，重复2，3步。
5、反复重复1-4步，得到一系列弯曲射线路径和相应的旅行时序列T0＞＝T1＞＝…，如果某个旅行时TM满足TM-1-TM＜＝ε，则追踪停止。
该技术特点是，能避免射线盲区和射线追踪路径并非时间最短路径问题，但此方法存在两点技术问题1、射线追踪速度慢在具体操作中，在作扰动的过程中，必须对所有可能的点依次进行扰动，这样就增加了运算量，大大降低了射线追踪的速度。
2、扰动量的选择问题在实际的操作和应用中，若扰动量选择过�。纫话憧梢源锏揭螅氖泵飨栽黾�。反之，若扰动量选取过大，则达不到最小旅行时的要求，因此需要对扰动量进行反复多次试验，以确定有效和合理的取值范围。以达到一方面即可以满足计算精度的要求，另一方面又不增加无谓的运算量。在实际的应用中，扰动量的合理选取十分困难。
发明内容
现有技术中的矩形网格三点扰动射线追踪法是针对速度模型离散的情况提出的，不论什么样的速度场、低速带的影响等，总能找到最短时间路径，避免了射线盲区和追踪路径并非时间最短路径问题，是一个较为理想的射线追踪方法。但该方法存在两大难题，其一是速度太慢。在作扰动的过程中，必须对所有可能的点依次进行扰动。然后计算旅行时，进行最小旅行时的比较，确定最终的路径的位置。其二是扰动量大小的取舍问题。扰动量过大，速度能加快，但找不到最短时间路径；扰动量过�。扑憔却锏胶螅钚÷眯惺辈换嵊腥魏问导室庖宓母谋洌俣让飨约趼�。只有一个合适的扰动量才能即满足计算精度，又不增加无谓的运算量。最令人困扰的还有这个扰动量随计算深度的不同，需要进行调整。深度增加时，扰动量需要适当地减小。在实际计算过程中，实现这种要求相当困难。
本发明是基于矩形网格三点扰动射线追踪法，利用Fermat最小旅行时原理，直接取得扰动中间点的位置，不需要对中间点进行依次扰动，从而提高了矩形网格三点扰动射线追踪法计算速度，同时由于在Fermat最小旅行时原则下计算最小旅行时的位置，因此避免了扰动量选取问题。
具体的方法如下一种用于地震勘探中射线追踪的方法，其特征在于本方法将射线追踪场进行矩阵网格化后，利用矩形网格三点和FERMAT原理，逐一对射线路径点进行追踪，形成最短走时射线追踪路径；其具体方法为A、对射线追踪区域进行网格化操作；B、建立深度域的二维速度模型的离散速度�。籆、对离散速度场做平滑处理每个深度界面的网格，沿水平轴方向取奇数个点，如5，7，9等进行数学平均或中值滤波的平滑处理；D、设定初始射线路径和该初始路径与所述网格的交点p1、p2、…pN::
E、计算初始射线路径的旅行时间T0；F、逐点追踪射线路径(1)取点取初始射线路径上三个顺序交点(p1，p2，p3，)；(2)计算固定两端点p1和p3，对中间点p2，用Fermat原理计算得到最小旅行时对应的中间点q1，并用新中间点q1代替原有的中间点p2，使p1q1p3为p1到p3的时间最短的弯曲射线路径；(3)再取q1，p3，p4三点，对其进行修正，重复(2)的操作，这样依次取三点修正直到修正到射线路径的终点pN为止；(4)经过步骤(1)-(3)后，形成一条新的射线路径p1q1q2…pN，(5)计算旅行时间根据p1q1q2…pN路径计算新路径的旅行时间T1；(6)阈值比较将前后两次旅行时间(T0和T1)的差值和阈值进行比较；(7)重复步骤若阈值比较的结果不符合精度要求，则以p1q1q2…pN作为初始路径，重复上述步骤F中的(1)-(6)的追踪过程，得出新旅行时间T2，并与T1进行阈值比较，不满足精度要求则反复重复步骤F中的(1)-(6)的追踪过程，不断形成新的射线路径，逐点逼近最短走时的射线路径；G、输出射线路径步骤(1)若步骤F中的(6)阈值比较的结果符合精度要求，即旅行时间的差值符合精度要求|TM-1-TM|＜＝ε，追踪停止；(2)输出相应的射线路径即为最小旅行时间的地震射线追踪路径，输出此射线路径。
在具体的应用中，所述步骤A的网格化操作中横向网格间隔大于纵向网格间隔。
在具体的应用中，所述步骤A的网格化操作中网格间隔为横向网格间隔与纵向网格间隔的比例为1∶1。
在具体的操作中，所述步骤B中的建立深度域的二维速度模型的离散速度场为从地震资料中，用速度分析方法得到纵波和对应深度的横波的叠加速度。用狄克斯公式将叠加速度换算为层速度，并将时间-速度域转换到深度-速度域，形成深度域的速度模型。在震源点-检波点之间对速度场进行内插和外插就生成深度域射线追踪需要的二维速度模型的离散速度场。
在具体的操作中，所述步骤C中的建立深度域的二维速度模型的离散速度场操作过程中作平滑处理，取5-9道的平滑。
在具体的应用中，所述步骤F中的FERMAT原理，其核心公式t(r)=(x-r)2+h12v1(r)+r2+h22v2(r)dt(r)dr=r-xv1(x-r)2+h12+rv2r2+h22=0]]>其中x网格交点Pi+2-Pi水平间距；r网格交点Pi+2-Pi+1水平间距；h1网格交点Pi+1-Pi垂直间距；h2网格交点Pi+2-Pi+1垂直间距；v1，v2分别为Pi+1点上下网格速度；t从Pi～Pi+2的旅行时；dt/dr时间对r的求导。
在实际的应用中，所述步骤F中的逐点追踪递推公式x＝p3-p1、q1＝p1+(x-r)。
x＝pi+2-qi、qi+1＝qi+(x-r)，(i＝1，2，……，N-1)其中r为Fermat原理求解的具有最小旅行时的位置。
式中p1、p2，…pN系列为初始网格交点的水平坐标；q1是对p2点修改后的水平坐标；然后对p系列依次进行修改，并利用前面修改后的q值，形成新的q系列，也就得到了p1至pN新的射线路径。
根据上述的技术方案，本申请的矩形网格三点Fermat射线追踪法是基于矩形网格三点扰动法的一种提高计算速度的方法。取矩形网格三个点，在Fermat最小旅行时原则下求取中间点的位置，而不象扰动法那样依次扰动。因此，计算速度比扰动法提高2倍多，同时不受扰动过程中的扰动量大小选择的困扰。该方法继承了矩形网格三点扰动法的优点，对任意离散的速度�。苣苷业阶疃淌奔渎肪叮苊饬松湎呙で妥纷俾肪恫⒎鞘奔渥疃搪肪段侍�。在地震偏移、反演和成析成像方法应用中具有广泛的应用前景。
图1为本发明从S点到R点，穿过界面的射线旅行示意图；图2为本发明矩形网格三点Fermat射线追踪具体算法示意图；图3为本发明矩形网格三点Fermat射线追踪方法的流程图；
图4为本发明中P波深度域的层速度模型；图5为本发明50×50的网格三点Fermat射线追踪示意图；图6为本发明矩形网格三点Fermat射线追踪法的PP波最小旅行时射线路径图；图7为本发明和矩形网格三点扰动法射线追踪运行现场的操作界面示意图和耗时对比图；图8为本发明和矩形网格三点扰动法射线追踪方法的最小旅行时路径叠合显示图。
具体的实施方式图1为本发明从S点到R点，穿过界面的射线旅行示意图。在矩形网格三点扰动追踪法的基础上，本发明提出矩形网格三点Fermat射线追踪法，同样取矩形网格和射线路径的三个交点xi-1，xi和xi+1，对中间点xi，不是给出一个扰动量对其扰动，而是根据Fermat最小旅行时原理，确定具有最小旅行时的位置(r)。这样一次就能够确定出最小旅行时的位置。不是最小旅行的中间点，不必依次扰动。以新的中间点位置替代xi，顺序取下一个交点，反复上述过程，直到最后一个交点，得到新的射线路径，计算新路径的旅行时，以新的射线路径作为新的初始射线路径，重复上述过程，直到旅行时的改变量小于给定的精度，终止追踪计算。其计算速度上比扰动法快，而且继承了扰动法的优点，适应各种复杂速度场。
t(r)=(x-r)2+h12v1(r)+r2+h22v2(r)dt(r)dr=0]]>其中，x网格交点Pi+2-Pi水平间距；r网格交点Pi+2-Pi+1水平间距；h1网格交点Pi+1-Pi垂直间距；h2网格交点Pi+2-Pi+1垂直间距；v1，v2分别为Pi+1点上下网格速度；t从Pi～Pi+2的旅行时；dt/dr时间对r的求导。
如图1所示，上式中x＝xi+1-xi-1是固定两端点S(xi-1)和R(xi+1)的水平距离；r＝xi+1-xi是中间点(xi)到R(xi+1)的水平距离；v1、v2分别是扰动点上下网格的速度；h1，h2是扰动点上下网格的高度。t(r)是从S经过中间点xi，到达R的射线旅行时间，如果希望从S到R这段旅行时最�。匦氯范ㄖ屑涞鉿i的位置，三点扰动法是用一个扰动量在界面上不断改变xi位置，计算相应的旅行时t(r)，比较各xi的旅行时，求出最小旅行时，从而确定xi的位置作为最小旅行时的射线路径。发明的三点Fermat方法，不是对xi进行依次扰动，而是对旅行时t(r)求关于r的导数，当导数为零时，求解出关于r的解，对应的xi就是最小旅行时的射线路径。这就是Fermat原理。
图2为本发明矩形网格三点Fermat射线追踪具体算法示意图。
具体方法图中，a为初始射线，b为最小时射线。1、连接S，R的直线作为初始射线路径，求出此直线与网格线的交点，P1，P2，…，PN，并计算旅行时T0。
2、用P1，P2，P3三点修正初始射线路径，其中P1，P3保持不动，用Fermat最小旅行时公式，直接计算得出Q1点的位置，使P1Q1P3为P1到P3的时间最短的弯曲射线路径。
3、取Q1，P3，P4三点作修正，同第2步，这样依次取三点修正直到修正到PN为止，并计算新路径P1Q1Q2…PN的旅行时T1。
4、把PiQ1Q2…PN的连线当作新的初始射线路径，重复2，3步。
5、反复重复1-4步，得到一系列弯曲射线路径和相应的旅行时序列T0＞＝T1＞＝…，如果某个旅行时TM满足TM-1TM＜＝ε，则追踪停止。
图3为本发明矩形网格三点Fermat射线追踪方法的流程图，其具体步骤同上。
图4为本发明中P波深度域的层速度模型图。
取射线追踪的区域为深度0-5000米，最大震源点-检波点距离4000米。按照图3的计算流程，1、用40×50米的网格对射线追踪的区域进行网格化；2、从地震资料中，用速度分析方法可以得到纵波(P波)和横波(S波)的叠加速度。用狄克斯公式将叠加速度换算为层速度，并将时间-速度域转换到深度-速度域，形成深度域的速度模型。在震源点-检波点(接收点)之间对速度场进行内插和外插就生成深度域射线追踪需要的二维速度模型的离散速度场(见图4)。离散速度场的速度在纵、横向变化明显，速度层中夹有高速层，高速层下有低速度层，整个速度场比较复杂；3、给出初始射线路径，计算初始射线路径与网格交点，并计算初始射线路径的旅行时；4、循环交点，取射线上三个顺序的交点，固定两端点，对中间点用Fermat原理计算最小旅行时对应的中间点，用新的中间点代替老的中间点；5、判断交点是否到最后一个点？否，用新的点和顺序向下取射线交点，组成新的三个点，重复步骤步3。是，计算新路径的旅行时；6、判断旅行时是否满足计算精度？否，以新的路径当作新的初始射线路径，重复步骤3、4和5。是，终止射线追踪。
图5为本发明50×50的网格三点Fermat射线追踪示意图；用50×50的网格，网格间距(50米)，速度场在网格作横向和纵向线性变化。固定两点(震源点-检波点)，进行初始射线追踪和最小旅行时追踪。图5是震源点-检波点的两点射线追踪的射线图。直线是初始射线，曲线是最小旅行时射线。弯曲射线的旅行时最�。鄙湎卟⒎鞘亲钚÷眯惺钡穆肪�。由于速度场线性变化，所以弯曲射线十分圆滑。其中，a为初始射线，b为最小时射线。
图6为本发明矩形网格三点Fermat射线追踪法的PP波最小旅行时射线路径图。从深度点2000米至4800米，固定震源点-反射点-检波点，矩形三点Fermat射线追踪PP波射线路径图。从图6可见，射线除了明显弯曲外，没有出现盲区，所有射线都能够通过复杂速度模型，以最小旅行时到达地面接收点。
图7为本发明和矩形网格三点扰动法射线追踪运行现场的操作界面。矩形网格三点扰动法射线追踪运行现场(上)，输入扰动量为0.05，运行开始时间为8点58分41秒，深度点为2500米，PP波(PP波为纵波)的最小旅行时为1814毫秒，PS波(PS波为转换波)的最小旅行时为2522毫秒，终止时间为8点59分24秒，计算耗时43秒；本发明矩形网格三点Fermat射线追踪法运行现场(下)，运行开始时间为9点4分9秒，深度点为2500米，PP波的最小旅行时为1812毫秒，PS波的最小旅行时为2520毫秒，终止时间为9点4分29秒，计算耗时20秒。
在SGI工作站上，用矩形网格(网格间隔40×50)三点扰动法，扰动量选0.05(扰动量小于0.05时，最小旅行时的改变量不到4毫秒，但耗时明显增加。扰动量过大，如取0.1，则达不到最小旅行时的要求)，因此，扰动量0.05是经过试验后确定的量，可以即满足计算精度，又不增加无谓的运算量。在震源点，P波从三个角度入射到地下界面，在深度点2500米的界面反射P波和转换S波到地面检波点，PP波最小旅行时为1814毫秒，PS波最小旅行时为2522毫秒，计算总耗时为43秒，运算现场如图7上。同样的计算精度(甚至比扰动法精确2毫秒，图7下)、同样的深度点和震源点-检波点，矩形网格三点Fermat射线追踪法的计算总耗时为20秒。提高速度大约2.15倍。两种方法的PP波和PS波最小旅行时的射线路径见图8，实线是矩形网格三点扰动法的纵波射线追踪的结果，空心圆为矩形网格三点Fermat法纵波射线追踪法的结果；小方框是矩形网格三点扰动法的转换波射线追踪的结果，十子为矩形网格三点Fermat法转换波射线追踪法的结果。PP和PS波的下行入射P波路径相同，在界面上遵守Snall定律，反射P波和转换S波上行至检波点位置。两种方法的射线路径基本重叠。从最小旅行时的结果也可以看出，2毫秒的差别在射线路径上基本反映不出多少差别。但两种方法的运算速度差别很大，在深度2500米处，两者的计算速度差2.15倍，在浅层，如1000米处，计算结果表明，两者的计算速度差2.67倍。即在相同的扰动量条件下，界面愈浅，两种方法的计算速度差别愈大。
对于矩形网格三点扰动算法，随着深度的增加，固定一个扰动量还不能满足最小旅行时的要求。需要根据深度的增加，适当减小扰动量，这样使这种方法的耗时进一步增加。因此矩形网格三点Fermat射线追踪法比三点扰动法运算速度提高2倍多，在实际的应用中得到广泛的应用。
图8为两种射线追踪方法的最小旅行时路径叠合显示图。
图中pp-r三点扰动法的纵波射线；pp-f三点Fermat法的纵波射线；ps-r三点扰动法的转换波射线；ps-f三点Fermat法的转换波射线。
权利要求
1.一种用于地震勘探中射线追踪的方法，其特征在于所述方法包括以下步骤A、对射线追踪区域进行网格化操作；B、建立深度域的二维速度模型的离散速度�。籆、对离散速度场做平滑处理对于每个深度界面的网格，沿水平轴方向取奇数个点进行数学平均或中值滤波的平滑处理；D、设定初始射线路径和确定该初始路径与所述网格的交点p1、p2、…pN；E、计算初始射线路径的旅行时间T0；F、逐点追踪射线路径(1)取点取初始射线路径上三个顺序交点(p1，p2，p3，)；(2)计算固定两端点p1和p3，对中间点p2，用Fermat原理计算得到最小旅行时对应的中间点q1，并用新中间点q1代替原有的中间点p2，使p1q1p3成为p1到p3时间最短的弯曲射线路径；(3)再取q1，p3，p4三点，对p3进行修正，重复(2)的操作，这样依次取三点修正直到修正到射线路径的终点pN为止；(4)经过步骤(1)-(3)后，形成一条新的射线路径p1q1q2…pN；(5)计算旅行时间根据p1q1q2…pN路径计算新路径的旅行时间T1；(6)阈值比较将前后两次旅行时间(T0和T1)的差值和阈值进行比较；(7)重复步骤若阈值比较的结果不符合精度要求，则以p1q1q2…pN作为初始路径，重复上述步骤F中的(1)-(6)的追踪过程，得出新旅行时间T2，并与T1进行阈值比较，不满足精度要求则反复重复步骤F中的(1)-(6)的追踪过程，不断形成新的射线路径，逐点逼近最短走时的射线路径；G、输出射线路径步骤(1)若步骤F中的(6)阈值比较的结果符合精度要求，即旅行时间的差值符合精度要求|TM-1-TM|＜＝ε，追踪停止；(2)输出相应的射线路径即为最小旅行时间的地震射线追踪路径，输出此射线路径。
2.根据权利要求1所述的方法，其特征在于所述步骤A的网格化操作中横向网格间隔大于纵向网格间隔。
3.根据权利要求1或2所述的方法，其特征在于所述步骤A的网格化操作中横向网格间隔与纵向网格间隔的比例为1∶1。
4.根据权利要求1所述的方法，其特征在于所述步骤B中的建立深度域的二维速度模型的离散速度场从地震资料中，用速度分析方法得到纵波和对应深度的横波的叠加速度；再用狄克斯公式将叠加速度换算为层速度，并将时间-速度域转换到深度-速度域，形成深度域的速度模型；然后在震源点-检波点之间对速度场进行内插和外插就生成深度域射线追踪需要的二维速度模型的离散速度场。
5.根据权利要求1所述的方法，其特征在于所述步骤C对离散速度场作平滑处理中取5-9道的平滑。
6.根据权利要求1所述的方法，其特征在于所述步骤F中的Fermat原理，其应用公式为t(r)=(x-r)2+h12v1(r)+r2+h22v2(r)]]>dt(r)dr=r-xv1(x-r)2+h12+rv2r2+h22=0]]>其中，x网格交点Pi+2-Pi水平间距；r网格交点Pi+2-Pi+1水平间距；h1网格交点Pi+1-Pi垂直间距；h2网格交点Pi+2-Pi+1垂直间距；v1，v2分别为Pi+1点上下网格速度；t从Pi～Pi+2的旅行时；dt/dr时间对r的求导。
7.根据权利要求1所述的方法，其特征在于所述步骤F中的逐点追踪递推公式x＝p3-p1、q1＝p1+(x-r)；x＝pi+2-qi、qi+1＝qi+(x-r)，(i＝1，2，……，N-1)；其中r为Fermat原理求解的具有最小旅行时的位置；p1、p2，…pN系列为初始网格交点的水平坐标；q1是对p2点修改后的水平坐标；然后对p系列依次进行修改，并利用前面修改后的q值，形成新的q系列，也就得到了从p1至pN新的射线路径。
全文摘要
本发明属于石油地震勘探资料处理方法技术领域，尤其涉及一种最小旅行时地震射线追踪的方法。本发明的矩形网格三点Fermat射线追踪法是基于矩形网格三点扰动法的一种提高计算速度的方法。取矩形网格三个点，在Fermat最小旅行时原则下求取中间点的位置，而不象扰动法那样依次扰动。因此，计算速度比扰动法提高2倍多，同时不受扰动量大小选择的困扰。该方法继承了矩形网格三点扰动法的优点，对任意离散的速度场，总能找到最短时间路径，避免了射线盲区和追踪路径并非时间最短路径问题。在地震偏移、反演和成析成像方法应用中具有广泛的应用前景。
文档编号G01V1/28GK1712991SQ200410049719
公开日2005年12月28日 申请日期2004年6月25日 优先权日2004年6月25日
发明者黄中玉 申请人:中国石油化工股份有限公司, 中国石油化工股份有限公司石油勘探开发研究院
专业数控铣床产品供应商
20年生产经验，实力工厂，精工品质，质量有保障
