专利名称：：一种分形地震波初至拾取的方法
技术领域：
：本发明涉及一种在油气勘探中自动拾取地震初至波的技术，具体地说涉及一种分形地震波初至拾取的方法。
背景技术：
：自从Mandelbrot(1967)在Sciences上发表了名为英国的海岸线有多长的著名文章之后，分形理论开始得到广泛地研究。一个几何形体的维数不再局限于整数，而是扩展到了分数范围。Mandelbort(1982)在他的经典著作中描述了海岸线的测量方法之后，分形学以海岸线的几何形态为例，发展了关于曲线的分维的计算方法(FalconerK.J.,1990)。近年来，分形学渗透到地震勘探研究领域，地震记录可以被看作为一条曲线，利用改变粗视化(roughness)程度的方法(Feder，1988)可以对地震记录进行分维计算，根据分维的变化特征可以判别地震波的初至走时(Boschettietal.,1996)。地震记录的维数被定义在大于1小于2的范围内。利用地震道分维的变化判别地震初至波的依据是，地震记录上初至波到达之前，分维数值为常数，从初至波到达时刻起，分维数值开始发生变化(Boschettietal.,1996)。但是当地震记录存在一定的噪声时，地震记录的分维在初至波到达之前不可能是常数。因此，分形算法在遇到强噪声背景地震记录时，存在如何提高算法的抗噪声能力的问题。本发明针对这一问题引进广义分维的概念，大大提高了分形算法的抗噪声能力。在经典分形理论中，地震记录的分维变化区间是[1，2]，在广义分维的范畴，地震记录分维数值的变化区间可以扩展到[O,>2)。此时维数不再有物理意义，而仅仅指示了地震记录在某种测量方式下的相对特征值。这一相对特征值为本发明拾取高噪声背景下的地震初至波走时，甚至反射波走时，提供了极为方便的定量依据。准确拾取地震波走时是地震勘探以及地震层析成像方法中计算速度的重要依据。地球表面的未压实层和风化层相对于较深部的固结层和岩石层的速度存在相当大的差异，表层厚度及速度的估算直接影响其下覆地层速度估算的准确性，所以正确拾取地震波初至走时在表层厚度及速度计算中意义重大。由于表层地质结构的复杂性及人们对表层估算精度要求的不断提高，对地震波初至的理解已不再局限于直达波，初至是直达波与折射波的共存与复合已成为人们的共识。初至走时的计算机自动拾取方法也因此而不断改进。近年来，许多文献记载了初至走时自动拾取的方法。例如Gelchinsky等用相邻道互相关方法，Coppens提出在不同大小的时窗内进行能量比较的方法Ramananantoandro用改变褶积算子宽度直到振幅最大值的位置不再发生变化的方法，Boschetti等于1996年发表了基于分形理论的dividermethod和Hurstmethod可用来检测地震波的初至走时的文章等。大多数初至走时自动拾取方法的关键是依据初至与其相邻点之间振幅及频率的差异。
发明内容本发明的目的在于提供一种分形地震波初至拾取的方法。本发明提供了地震道Hausdorff广义分数维的算法以及如何根据分数维的变化标定地震波初至走时的方法。地震波初至到达之前，地震道时间序列为随机噪声的响应，地震波初至到达之后，地震道时间序列为随机噪声与地震波动的复合。本发明在分形学理论的基础上利用Hausdorff广义分数维识别地震道初至走时不仅是一种快速、准确的方法，而且在计算中不需人工干预，能够实现计算机全自动拾取。这种方法对相邻道无特殊限定，适用于各种不规则的二维及三维地震勘探观测系统、人工地震透射波层析成像的走时识别，也适用于天然地震层析成像研究中地震波初至走时的识别。本发明提供的分形地震波初至拾取的方法，步骤如下a)首先对地震记录进行去噪处理，消除环境噪声和有规律的系统噪声；b)定义一个工作窗，工作窗的长度要保证有足够的样本量计算分形维数，以及计算速度；c)选取一个任意半径r的圆作为地震记录度量的初始尺度，以地震记录上第一个采样点为圆心，沿地震记录移动圆；d)通过二维圆方程与地震曲线方程的联合求解，计算圆与地震记录曲线的交点；e)以该交点作为圆度量地震曲线长度的下一个起始点，同样以r为半径画圆，重复步骤c的操作，再次计算出圆和地震记录新的交点，直到将圆移动到地震记录的末端，将地震记录长度计算出来；f)给圆半径r一个增量，以新的圆半径重复步骤b、c、d的操作，计算出一个新的地震记录长度丄；g)当一个工作窗的曲线段完全被圆覆盖后，根据公式D^"A"，'计算该工作窗的广义分形维数；h)移动一个数据采样点的距离，开始新一个工作窗的分维计算。i)重复步骤d、e、f、g、h操作，得到整个地震记录的广义分维数值；j)根据分维数值由常量变为变量的位置确定地震波初至，并将这一位置对应的地震波传播时间写入初至走时数据文件。本发明可用于地震资料的爆炸震源、锤击震源、伪随机可控震源和微震震源的初至波的自动拾取。拾取过程不需要人工干预，适用于山区崎岖地形的人工地震数据和非规则观测系统的人工地震数据以及天然地震观测数据。本发明相对公知技术而言所具有的优点和效果本发明提出地震初至波的自动拾取广义分维方法，广义分维方法具有定量化程度高、抗噪声能力强，并且不受观测系统约束，可以在单个地震道上实现的优点。这种方法不仅适用于规则观测的人工地震记录，而且适用于非规则观测的天然地震记录。广义分维的引进提高了分形算法初至波自动拾取的抗噪声能力，为解决低信噪比地震记录的走时拾取问题提供了一种实用方法。对于人工地震勘探中可控震源的地震记录，可以首先将可控震源地震记录作最小相位处理，使地震信号的能量集中在子波的最前部，使之与炸药震源的信号相似，即可望得到与炸药震源地震记录同样的拾取效果。图1是地震道的分维计算.图2(a)是地震道，长方形表示工作窗；图2(b)是拟和计算的分维曲线。图3是圆的半径与工作窗长度关系;如果圆的半径大于工作窗的长度，计算的地震记录长度近似为零。图4(a)-(c)是不同小圆半径计算的地震记录长度不同。图5(a)-(i)是不同测量尺度下，图中(a)所示地震记录分维数值发生的变化。图6(a)-(e)是在地震记录(a)上加入不同等级噪声，观察分维算法的抗噪声能力。(a)表示一个实际的地震记录和它的广义分维曲线。(b)-(e)分别表示增加了20%、40%、60%、80。/。不同程度随机噪声的地震记录及其广义分维曲线。在(a)-(e)中，每一条靠下面的曲线是地震记录，上面的曲线是其广义分维曲线。A表示振幅，D表示维数，A和D为无量纲数。图7(a)-(b)是对中国华北遥测地震台网观测到的地震波形记录用广义分维算法获得的P波震相走时的自动拾取结果。(a)和(b)分别表示两个地震记录及其广义分维曲线，靠下的曲线是地震记录，靠上的曲线是其分维曲线。图中A表示振幅，D表示分维数值，A和D是无量纲数。图8(a)-(b)是对一个复杂地形地震探区的三维地震数据用广义分形算法获得的地震P波初至走时的自动拾取结果。即，是对一个复杂地形地震探区的三维地震数据用广义分维获得的地震P波初至走时的自动拾取结果。图9(a)-(d)是从图8中选择的不同噪声背景的地震记录及其广义分维曲线。具体实施例方式地震记录是由噪声和信号共组成的，当地震波初至到达之前，地震记录为噪声背景，在初至到达之后，地震记录是噪声与信号的叠加。可以设想噪声是一个"物体"，而噪声加信号是另外一个"物体"，这两个"物体"有不同的几何特征，因此对这两个"物体"同时进行观测时，在广义分维范畴，可以选择一种合适的尺度，使对噪声维数的观测结果为零，对信号加噪声部分的维数的观测结果不为零。这个合适的尺度的选择可以通过试验得到。一旦这个合适的尺度(临界值)找到了，就可以用广义分维的概念来监测地震波初至的位置。(1)首先对地震记录进行去噪处理，消除环境噪声和有规律的系统噪声。(2)定义一个工作窗的长度如图1，工作窗的长度既要保证有足够的样本量计算分形维数，又要考虑计算速度。(3)选取一个半径为r的小圆作为地震记录度量的初始尺度，以地震记录上第一个釆样点为圆心，度量尺度r为小圆的半径，沿地震记录移动小圆。(4)通过二维圆方程与地震曲线方程的联合求解，计算圆与地震记录曲线的交点。(5)以新的交点作为小圆度量地震曲线长度的下一个起始点，以r为半径画圆，如步骤(3)移动小圆，再次计算出圆和地震记录新的交点，直到将小圆移动到地震记录的末端，将地震记录长度计算出来。(6)给小圆半径r一个增量，以新的半径作为小圆半径重复步骤(2)、(3)和(4)，计算出一个新的地震记录长度丄。(7)当一个工作窗的曲线段完全被小圆覆盖后，根据DK"g^""r'式计算该工作窗的广义分形维数。(8)移动一个数据釆样点的距离，开始新一个工作窗的分维计算。(9)重复步骤(4)到(8)，完成整条曲线(即整个地震记录)的广义分维计算。(10)根据分维数值由常量变为变量的位置确定地震波初至，并将这一位置对应的地震波传播时间写入初至走时数据文件。实施例1图5(a)是一个实际的地震记录。在0-1.2sed的记录段，地震记录上主要为噪声背景，1.2sed之后是地震信号与噪声背景的叠加，为简便起见，称1.2sed之后为信号段。对这个地震记录用不同的测量尺度按照图4所示方法计算它的分维曲线。图5的(b)、(c)、(d)、(e)、(f)、(g)、(h)、(i)给出了不同测量尺度下，图5(a)所示地震记录分维数值发生的变化。对图5的结果作如下分析。当工作窗的长度(采样点数)一定时，测量尺度r的取值范围决定了地震记录分维曲线的特征(l)r在每一个工作窗长度内的取值范围集中在较大值附近时，地震记录噪声段的分维数值低于信号段，如图5的(b)、(C)、(d)。(2)r在每一个工作窗长度内的取值范围集中在较小值附近时，地震记录噪声段的分维数值高于信号段，如图5的(g)、(h)、(i)。(3)r在每一个工作窗长度内的取值范围集中在中间值附近时，地震记录噪声段与信号段的分维数值无明显差异，如图5的(e)、(f)。图5中(b)和(c)所示分维曲线的最大值都超过了2，因此，图5(b)和(c)所用的测量尺度在广义分维概念下成立。此外，r在每一个工作窗长度内的取值范围越大，取值数量越多，拟合Mandelbrot-Richardson曲线的样本量越大，地震记录的分维曲线越光滑，反之越不光滑。但样本量的大小只改变分维的数值，不改变噪声段与信号段变化点的分辨能力。在地震记录初至震相的自动拾取技术中，并不需要利用分维数值的稳定性，而只是利用连续分维曲线发生突变的位置，因此在这个问题上，样本量的大小并不重要，弄清这一点是节省计算量的关键。实施例2为了考察广义分维对低信噪比地震记录的适应能力，选择了一个实际的地震记录，对这个地震记录分别加入了振幅最大值的20%、40%、60%和80%几种不同程度的随机噪声，图6给出了考察计算结果。由图6可以看出，广义分维算法具有极强的抗噪声能力，而经典分形算法则不具备这样强的抗噪声能力。图6表明，对不同信噪比的地震记录，分维的数值发生变化，但特征保持不变。计算机可以容易地识别分维曲线的这一特征，从而实现初至波和后续强地震信号的自动拾取。实施例3用广义分维算法，本发明对中国华北遥测地震台网观测到的地震波形记录进行了P波震相走时的自动拾取，图7表示其中2个记录及其广义分维曲线。从图7可以看出，与每一个地震记录对应的广义分维曲线在P波到达之前表现为常数，从p波到达的时刻起发生剧烈变化。广义分维曲线发生变化的位置对应p波初至。利用计算机可以容易地判断地震记录的广义分维由常量变为变量的位置，从而对初至信号实现计算机的自动识别。与常规的初至自动拾取所用的互相关算法不同，广义分维算法不依赖于相邻地震记录的相似性，而是在单个的地震记录上实现，这对于天然地震记录非常适用。从图7中还可以看出，每一条地震记录都有各自的噪声背景和波形特征，但只要广义分维的测量尺度选择得合适，这些地震记录的广义分维都在p波震相出现之前为常数。实施例4本发明对一个复杂地形地震探区的三维地震数据进行了初至波的自动拾取。图8表示对该区一个炮集的拾取结果。图8中黑实线表示用广义分维全自动算法拾取到的初至波的位置。在图8中选择出不同噪声背景的地震记录，把这些地震记录和它们的分维曲线表示在图9中。由图9可见，尽管地震记录表现出不同的噪声背景，但只要测量分维所用的尺度的最小值取得合适，就可以获得能被用来进行初至波自动拾取的分维曲线特征。以上实施例说明，不论地震记录的噪声程度如何，它们的分维曲线都表现出在初至波到达之前保持为常数，在初至波到达之后发生剧烈变化的特征。这些结果表示地震记录广义分维特征可以用来进行快速、自动的初至波拾取。实施例5表l显示了三条地震记录的分维数值和初至位置。表l:地震道振幅值与地震道分形维数值<table>tableseeoriginaldocumentpage12</column></row><table>权利要求1、一种分形地震波初至拾取的方法，步骤如下a)首先对地震记录进行去噪处理，消除环境噪声和有规律的系统噪声；b)定义一个工作窗，工作窗的长度要保证有足够的样本量计算分形维数，以及计算速度；c)选取一个任意半径r的圆作为地震记录度量的初始尺度，以地震记录上第一个采样点为圆心，沿地震记录移动圆；d)通过二维圆方程与地震曲线方程的联合求解，计算圆与地震记录曲线的交点；e)以该交点作为圆度量地震曲线长度的下一个起始点，同样以r为半径画圆，重复步骤c的操作，再次计算出圆和地震记录新的交点，直到将圆移动到地震记录的末端，将地震记录长度计算出来；f)给圆半径r一个增量，以新的圆半径重复步骤b、c、d的操作，计算出一个新的地震记录长度L；g)当一个工作窗的曲线段完全被圆覆盖后，根据公式D∝LogLi/Logi计算该工作窗的广义分形维数；h)移动一个数据采样点的距离，开始新一个工作窗的分维计算。i)重复步骤d、e、f、g、h操作，得到整个地震记录的广义分维数值；j)根据分维数值由常量变为变量的位置确定地震波初至，并将这一位置对应的地震波传播时间写入初至走时数据文件。全文摘要一种分形地震波初至拾取的方法，涉及一种地震道Hausdorff广义分数维的算法以及如何根据分数维的变化标定地震波初至走时的方法。地震波初至到达之前，地震道时间序列为随机噪声的响应，地震波初至到达之后，地震道时间序列为随机噪声与地震波动的复合。本发明在分形学理论的基础上利用Hausdorff广义分数维识别地震道初至走时不仅是一种快速、准确的方法，而且在计算中不需人工干预，能够实现计算机全自动拾取。这种方法对相邻道无特殊限定，适用于各种不规则的二维及三维地震勘探观测系统、人工地震透射波层析成像的走时识别，也适用于天然地震层析成像研究中地震波初至走时的识别。文档编号G01V1/28GK101452082SQ20071017879公开日2009年6月10日申请日期2007年12月5日优先权日2007年12月5日发明者刘伊克,旭常申请人:中国科学院地质与地球物理研究所