专利名称：一种三维地震偏移成像的定点分解方法
技术领域：
本发明属于勘探地球物理领域中地震偏移成像方面所研究的课题，地震勘探数据是通过现场仪器设备录入到磁性介质并且在地震数据处理软件平台中形成图像，本发明提出一种简便算法对地震反射波数据进行偏移成像处理计算，为了解地下介质结构从而准确地确定勘探井位或其它地下资源所处位置的切实可行的可靠方法。
背景技术：
利用地震反射波数据资料对地层结构进行成像是地震勘探的核心技术之一。这种称为偏移成像(migration imaging)的技术，很像医学中的CT成像，但是其复杂程度却是难以类比。地震偏移成像的困难在于资料数据量特别大；信号源不确定；介质结构复杂；地震波传播中干扰噪声强大，绕射、屏蔽、吸收等无法掌握的不利因素都对最终成像精度产生坏的影响。与此同时，就成本而言，一口勘探井的工业投入一般达到数千万圆人民币。因此工业上总是要求快速准确的偏移成像成果，作为部署勘探方案、确定井位的依据。因此三维地震偏移成像技术的进步是期待已久的，在地球物理学领域，三维地震偏移成像技术从理论到工业应用是长期以来受到普遍关注、吸引了大批优秀科学家和工程技术人员来从事这方面的研究。
近二、三十年来，借助于计算机技术的迅猛发展，在专门用于地震资料处理大型计算机系统或网络机群上已经建立了不少各具特色的软件平台，这是能源工业勘探中必不可少的基本的常规信号处理工具。但是，任何一种平台都迫切需要提高偏移成像处理的速度和精度，目前还没有一种平台能提供令人满意的三维地震偏移成像技术，大批的技术力量在致力于研究开发新的方法。当前，受到普遍重视的是三维隐格式有限差分偏移方法，这种方法的优点是适应地下速度的横向变化，精度高；缺点是计算量太大。
1、数据量问题对于10km×10km的叠后数据资料，如果线距和CDP间距为25m，采样间隔为4.00ms，那么在每个时间点上的数据达到400×400＝160000个，在这个时间点上需要求解的方程阶数为160000×160000，对于记录长度为6秒的资料需要求解上述阶数的方程3000-6000次。
2、计算效率问题由于数据量大，计算繁杂通常已知的求解方法无法胜任。现有的计算数学理论，一般只研究一次或几次求解的算法。对于一块很小面积的三维地震数据如果数千次地求解阶数如此巨大的方程，已经不再是效率和成本的问题了，是有没有工业应用的前景问题了。
3、方向畸变问题由于波场外推过程中数千次求解二阶波动方程，每次都使用已算出的数据，如果处理中不能控制好畸变，那么误差的传递积累，足以使得最终成果毫无实际意义。
所以长期以来，人们期待的是快速准确的计算技术，如果计算技术方面的问题解决了，成像方面的问题也就会变的常简单。现有已知的计算技术中分裂求逆是比较常用的，大致有交替方向隐式(ADI)法也称为沿X-Y方向分裂(Brown，1983等)、基于螺旋边界的降维LU分解(Rickett andClaerbout，1998)。后者是一种非线性方法，计算中无法避免混沌现象，边界处理也有问题。前者在外推过程中由于圆对称性被破坏，无法克服方向畸变使最终成果被歪曲。在后面的图例和应用例，将给出具体说明。
发明内容
1、发明目的本发明采用一种新的分裂逼近方法，实现了一种不完全Cholesky分解(因为完全Cholesky分解的系数矩阵是稠密的不能用于计算)，规定极少数位置上的元素非零，采用误差最小的原则计算出这些非零元，在保证系数矩阵高度稀疏条件下完成逼近，这与传统方法不规定非零元位置具有本质的差别。由于不采用沿方向分裂，避免了方向畸变，保证偏移成像的精度。
2、技术方案本发明是通过以下技术方案来实现的一种三维地震偏移成像的定点分解方法，其特征在于该方法主要按下列步骤进行a、首先将有限差分法偏移的波场外推方程的系数矩阵A做定点的块LU分解分解前矩阵为 其中Tk，Ek为子块矩阵，Tk为三对角矩阵，Ek为对角矩阵，作块两对角分解为A＝LLT+R；其中L为块下三角矩阵，并且在主对角线上有非零元，主对角线上非零矩阵的左侧除相邻列有非零元之外，其余位置为零矩阵；LT为L的转置矩阵；R是块不完全Cholesky分解的剩余矩阵；零矩阵不需要进行数值计算，因此计算仅涉及主对角线上的非零元和主对角线两侧的两个非元，这样得到系数矩阵A的定点LU分解或称定点BIC分解；b、然后对定点BIC分解的残差进行计算其剩余矩阵R＝RL+RU为一个在上、下三角部分各有沿对角线方向的一列非零元的对称矩阵，其中RL为上三角部分的非零元矩阵，Rf为下三角部分的非零元矩阵，用x＝L-TL-1b-Rx很容易得到一种迭代算法x(k+l)＝L-TL-1b-Rx(k)；c、根据实际资料可以确定主元优势，进行外推步长计算，选择合适的步长；d、在嵌入到地震数据处理软件平台中的软件模块中，使用上述步骤中得到的A的定点分解式、及所选定的步长，即可得到精度较高的图像。
在步骤a中R的元相对于主元是小量，并且可以控制。
在步骤c中如果希望收敛更快速，可以在采集和预处理阶段通过道内插，或选取合适的步长实现。
3、优点及效果本发明为解决长期困扰勘探地球物理工业界的三维地震偏移成像精度问题，采用了一种全新的分裂逼近技术——定点不完全Cholesky分解，实现了对三维地震数据的快速准确的波场外推处理。使三维有限差分法偏移成像技术实现了从理论研究到工业应用的转变。此前由于要么是计算量过于巨大以至于没法承受，要么是波场外推中产生不能接受的畸变，使有限差分法三维偏移长期停留在理论研究的状态。本发明在数据处理过程中完全保持了波场传播的圆对称性，使最终成像精度得到可靠保证，又由于计算效率比传统方法提高几十倍，从计算成本上保证了工业上的实用性。本发明为勘探地球物理工业解决了一个长期没法解决的技术难题。实际应用中已经能够证实，所采用的计算方法在辽河油田计算所资料处理生产中实施中具有明显的技术效果和显著的经济效益。
四
附图1为三维地震偏移方法的分类图；附图2为针对波场外推算子圆对称性畸变所作的各种校正图；附图3为CG软件偏移的脉冲响应图；附图4为波场外推圆对称性效果图；附图5为另一种波场外推圆对称性效果图；附图6为本发明的波场外推圆对称性效果图；附图7为Q271线偏移前叠加剖面图；附图8为Q271线偏移效果图；附图9为另一Q271线偏移效果图；附图10为本发明的Q271线偏移效果图；附图11为三种Q271线偏移效果对比图；附图12为0251线偏移效果图；附图13为另一0251线偏移效果图；附图14为本发明的0251线偏移效果图；附图15为三种0251线偏移效果对比图。
五具体实施例方式本发明主要是通过快速准确的计算技术，将地震反射波数据资料进行计算，得到的数据结果在地震数据处理软件平台中形成较为准确的图像，从而得到对地层结构的较为精确的勘探成果。现有的三维偏移成像方法有很多种，见附图1所示。但对于地震勘探形成地层结构的图像却很难得到满意的精度，主要是由于计算量过大、外推过程畸变误差的积累或者在计算中产生混沌现象，以及边界处理等方面的问题。目前较常采用的计算方法有交替方向隐式(ADI)法和螺旋边界非线性降维LU分解法，螺旋边界非线性降维LU分解法在计算中不能避免混沌现象，附图2为针对波场外推算子圆对称性畸变所作的各种校正，资料中采用9点、13点、17点、25点和六角形McCllalen变换。附图3为CG软件偏移的脉冲响应，从图中可以看出偏移的噪声严重。附图4为外国软件CG偏移的脉冲响应水平切片，从图中可以看出圆对称性明显畸变，偏移噪声严重。附图5为软件GR偏移的脉冲响应水平切片，从图中可以看出圆对称性明显畸变，偏移噪声严重。下面通过实施例对本发明做进一步说明实施例对于10km×10km的叠后数据资料，采用线距和CDP间距为25m，采样间隔为4.00ms，将得到的数据可以形成波场外推方程的系数矩阵如下 a、首先将波场外推方程的系数矩阵A做定点的LU分解，其中Tk，Ek为子块矩阵，Tk为三对角矩阵，Ek为对角矩阵作块两对角分解为A＝LLT+R；其中L为块下三角矩阵，并且只在主对角线上和主对角线左侧相邻一列有非零元，其余位置为零矩阵；LT为L的转置矩阵；R是块不完全Cholesky分解的剩余矩阵；零矩阵不需要数值计算，该矩阵计算仅涉及主对角线和主对角线左侧的两个非零元，于是得到波场外推方程的系数矩阵A的定点LU分解形式；b、然后对定点BIC分解的残差进行计算，本发明所做的定点BIC分解是一种特别稀疏的LU分解。其剩余矩阵R＝RL+RU为一个在上、下三角部分各有沿对角线方向的一列非零元矩阵，其中RL为上三角部分的非零元矩阵，RT为下三角部分的非零元矩阵，用x＝L-TL-1b-Rx很容易得到一种迭代算法
x(k+l)＝L-TL-1b-Rx9(k)；c、根据实际资料可以确定主元优势，进行外推步长计算，选择合适的步长；d、在嵌入到地震数据处理软件平台中的软件模块中，使用上述步骤中得到的A的定点分解式、及所选定的步长，即可得到精度较高的图像。
附图6所示的图像为模型数据，从图中能够清楚的看出其严格保持圆对称性，没有偏移噪声。
在步骤a中R的元相对于主元是小量，并且可以控制。
在步骤b中x(k+l)＝L-TL-1b-Rx(k)用以改进算法的精度。这是一种可以人为控制精度的方法，使用时可以在计算机资源与计算精度之间做权衡。这是现有的商业软件中没有的特点。
在步骤c中如果希望收敛更快速，可以在采集和预处理阶段通过道内插，或选取合适的步长达到预期目的。经过如上处理后，得到的数据在地震数据处理软件平台中形成了较为理想的偏移成果剖面(如附图10)。附图7为Q271线偏移前叠加剖面。附图8为CG软件Q271线偏移效果图，从图中可以看出圈出位置成像较差。附图9为GR软件Q271线偏移效果图，从图中可以看出圈出位置成像不理想。附图10为本发明定点分解法的Q271线偏移效果图，从图中可以看出圈出位置的图像成像清晰。附图11为三种Q271线偏移效果对比图，左边是GR软件成像效果图，中间为定点分解法成像效果图，右边是外国软件CG成像效果图，通过比较可以看出中间图像成像清晰，效果明显好于左、右两图像的成像效果。附图12为GR软件偏移结果，图像中方框圈出部位的小褶皱形态不符合地质规律，圆圈处连续性虽然还可以，但层次不够丰富。附图13为CG软件偏移结果，图像中方框圈出部位的右翼形态可靠，而左翼反射丟失，且圆圈处连续性差，横向难以追踪。附图14为本发明的定点分解法的偏移结果，图中可以看出圆圈内部位信噪比高，方框内部位成像准确。背斜形态可分辨，但连续性不太理想，圆圈处连续性分辨率明显改善，并且盖层断点干脆。附图15为三种0251线偏移效果对比图，左边是软件GR得到的图像，中间是定点分解法得到的图像，右边是外国软件CG形得到的图像，从图中可以看出矩形框内中间图像的对称性明显优于两边的图像，成像效果明显好于左、右两边的图像。
通过本技术方案的实施，得到的偏移成像成果能够较为准确的反映出地下介质的结构，为部署勘探方案、确定勘探井位和提供可靠的依据。
权利要求
1.一种三维地震偏移成像的定点分解方法，其特征在于该方法主要按下列步骤进行a、首先将有限差分法偏移的波场外推方程的系数矩阵做定点的块LU分解分解前矩阵为 其中Tk，EK为子块矩阵，Tk为三对角矩阵，Ek为对角矩阵，作块两对角分解为A＝LLT+R；其中L为块下三角矩阵，并且在主对角线上有非零元，主对角线上非零矩阵的左侧除相邻列有非零元之外，其余位置为零矩阵；LT为L的转置矩阵；R是块不完全Cholesky分解的剩余矩阵；零矩阵不需要进行数值计算，因此计算仅涉及主对角线上的非零元和主对角线两侧的两个非元，这样得到系数矩阵A的定点LU分解或称定点BIC分解；b、然后对定点BIC分解的残差进行计算其剩余矩阵R＝RL+RU为一个在上、下三角部分各有沿对角线方向的一列非零元的矩阵，其中RL为上三角部分的非零元矩阵，RT为下三角部分的非零元矩阵，用x＝L-TL-1b-Rx很容易得到一种迭代算法x(k+1)＝L-TL-1b-Rx(k)；c、根据实际资料可以确定主元优势，进行外推步长计算，选择合适的步长；d、在嵌入到地震数据处理软件平台中的软件模块中，使用上述步骤中得到的A的定点分解式、及所选定的步长，即可得到精度较高的图像。
2.根据权利要求1所述的一种三维地震偏移成像的定点分解方法，其特征在于在步骤a中R的元相对于主元是小量，并且是可以控制的。
3.根据权利要求1所述的一种三维地震偏移成像的定点分解方法，其特征在于在步骤c中如果希望收敛更快速，可以在采集和预处理阶段采取道内插，或选取合适的步长来实现。
全文摘要
本发明属于勘探地球物理领域中地震偏移成像方面所研究的课题，特别是涉及一种通过提出的简便算法将得到的地震反射波数据进行计算，然后嵌入到地震数据处理软件平台中形成图像，从而准确地确定油井或其它地下资源最佳位置的方法。本发明采用一种新的分裂逼近方法，实现了一种不完全Cholesky分解(因为完全Cholesky分解的系数矩阵是稠密的不能用于计算)，规定极少数位置上的元素非零，采用误差最小的原则计算出这些非零元，在保证系数矩阵高度稀疏条件下完成逼近。本发明主要解决了现有技术中计算量大、计算中无法避免混沌现象，边界处理不好等方面存在的问题，由于不采用沿方向分裂，避免了方向畸变，保证偏移成像的精度。
文档编号G01V1/28GK1888935SQ20061004721
公开日2007年1月3日 申请日期2006年7月17日 优先权日2006年7月17日
发明者石殿祥 申请人:石殿祥