槽波地震层析成像技术.docx

层析成像技术又称为CT技术。CT(Computerized Tomography)技术，早先是医学界用x射线或超声波为病人迅速诊断而发展起来的，是医学上的一种成像技术；后来，CT技术很快被引入地球物理学领域开始用以研究地球内部构造及天然地震等。英国Mason(1979)首先将这一技术用来研究Thoresby煤矿顶板的应力分布，随后在英国等国家得到越来越广泛应用。我国在80年代后期，西安煤炭研究分院、河北省煤炭研究所等单位开展了煤层槽波及坑道无线电波CT技术的研究与应用。1地震层析方法发展的回顾地震层析成像自引入到地球物理学领域以来,由于地震勘探中不同观测方式所具有的投影角度不同,解决问题的难易程度不同,计算机发展水平的限制以及理论基础(射线理论和波动理论)等方面的因素,近20多年来的发展主要遵循着以下三条主线 。成谷,马在田,耿建华,曹景忠地震层析成像发展回顾11 从透射层析发展到反射层析由于地震数据最通常的接收方式是地面接收,因此人们最先想到的是利用反射资料进行反射层析的研究。从80年代到现在为止,反射层析的研究一直没有中断过。例如,Cutler 和 Bishop的“地震层析成像:公式和方法”;Bube和 Resnick分析了反射地震层析成像能够测定的和不能准确测定的特征参数;Jannane等研究了从地震反射数据能够分辨的地球结构的尺度;Williamson针对反射层析中反射界面的情况提出了三种处理方式;Wang和Houseman利用反射地震振幅来反演速度的变化;Wang 和Pratt讨论了反射层析中地震旅行时和振幅对速度的不同分量的敏感度。与透射层析相比,反射层析成像除了受到“地震数据有限带宽、观测方式的有限视角”等问题的影响外,还受到“低速带、速度和反射界面深度之间的多解性”等特定因素的影响 ,而且反射层析与透射层析相比受有限视角的影响更大。因此,从总体上来说,层析成像的研究以透射层析(主要是跨井观测方式 ,也有 VSP 以及反向 VSP 观测方式)为多。在地震层析成像发展的早期 ,针对射线层析中大型 、欠定(或混定)、不相容、非常稀疏、有时病态或奇异的线性代数方程组而推出的各种解法,尽管对反射层析也同样适用 ,但基本上都是按照透射层析的情况进行试算的 。 针对层析成像中欠定(或混定)引起的解非唯一性问题而应用的各种平滑、约束和归一化的技巧 ,也主要是根据透射层析的情况进行研究的 ,如 Meyerholtz 和 Pavlis研究了地震层析中褶积压制的作用;Phillips和Fehler比较了各种平滑和归一化方式(阻尼最小平方法、阻尼最小平方加空间平均平滑器法、阻尼最小平方加均值平滑器法、褶积压制法等)应用于同一含噪音的合成数据的效果;Zhou在井间旅行时层析中采用了空间相干滤波器;Nemeth 和 Normark在井间旅行时层析中采用了动态平滑的措施。随着透射层析发展逐渐成熟 ,以及对反射层析中存在问题的相应解决办法的研究不断深入 ,人们发现 ,与透射数据相比 , 反射数据能够反映地下更精细尺度的结构和参数分布,并能节约观测成本,为此,反射层析的研究又重新提到日程上来 ,重新得到了人们的重视 。12 从旅行时层析发展到波形层析地震剖面主要包括旅行时信息和振幅信息。然而,同时利用以上两种信息对地下不同尺度的参数分布进行层析反演的结果很不理想 ,目标函数中局部极值的个数非常多,这意味着层析反演的非线性特强 。 因此人们转而开始研究地震剖面上的数据与地下介质模型的对应关系(抽象地说是数据空间与模型空间的对应关系)。J annan和 Beydoun等发现:数据空间与模型空间的对应关系并不简单。旅行时信息反映速度的大尺度变化 ,振幅信息反映速度的小尺度变化 ,速度的中尺度变化在数据中没有反映 ,但可以通过增大接收的排列长度和增加震 - 检系统的低频响应来缩小大尺度变化和小尺度变化之间的间距。受到这些研究结果的启发,人们将数据空间中的信息分解为旅行时和振幅 ,将模型空间中的慢度变化分解为长波长、中波长和短波长。在层析成像方法发展的初期 ,人们对地下介质参数变化的尺度分辨要求不高 ,仅利用旅行时信息进行层析反演已经足以反映地下介质参数变化的趋势,而且仅利用旅行时信息层析反演的结果还远较利用波场信息层析反演的结果为好。这是因为旅行时信息对应地下介质参数变化的趋势 ,因而目标函数中局部极值个数少,非线性程度相对减弱 ,使得计算和求解相对稳定。同时仅利用旅行时信息 ,意味着层析反演所用的数据量减少,计算量也相应减少。因此,尽管旅行时层析中旅行时拾取的精度对层析反演结果的影响非常大,而叠前剖面上信噪比低又给旅行时的拾取造成了困难 ,人们仍然对仅利用旅行时信息的层析方法情有独钟。然而,随着人们想要解决的地下介质情况越来越复杂,在叠前剖面上准确拾取旅行时越来越困难,这给旅行时层析的应用增添了一定的难度。同时人们也不再满足于旅行时信息所提供的地下介质参数分布的分辨率,而计算机的发展又使人们能够处理更多的数据量,因而利用地震剖面上更多的波场信息层析反演更精细尺度的地下介质参数的分布,已经成为发展的必然趋势。尽管如此,由于利用波场信息的局部寻优的层析方法对初始模型的依赖性强,而不依赖初始模型的全局寻优算法又计算量巨大,因而通常利用旅行时层析的结果作为波场层析的初始模型。13从射线层析发展到绕射层析在层析成像发展的早期,计算机水平较低,应用波动理论进行层析成像研究有一定的困难,因此人们转而利用波动理论的高频近似 射线理论 ,发展了很多节省计算机内存和运算量的方法。相对而言,射线层析利用的数据量少,计算效率高,计算相对稳定,因而受局部极值的影响较小。但由于受理论基础的限制,射线只能穿过物体的有限部分,对有限视角的敏感度高,同时由于利用的数据量少 ,层析反演出来的结果分辨率也就低,不能精细地解决地下介质的结构与参数分布情况。一般来说,应用射线层析要得到一个令人满意的结果要求在一个波长的范围内波速的变化很小,非均匀体的尺度与地震波长相比大很多,而且要求介质模型平滑。随着计算机能力的提高,基于波动理论的绕射层析方法发展了起来。绕射层析因为绕射作用,能够穿过射线层析中射线不能穿过的区域 ,使数据的覆盖面积扩大,但对有限视角的敏感度低。此外,绕射层析要求照明波长与非均匀体大小相当即可 ,对地下介质模型也没有特殊的要求。正是因为这些优点,地球物理绕射层析自 De2将声波绕射层析成像从常背景下推广到可变背景的情形;Beylkin和 Burridge将它扩展为弹性波绕射层析成像。针对以往的绕射层析成像法一般在空间2频率域进行 ,后来的一些研究者又阐述了空间域中的多频反投影重建问题;开发了层状背景的绕射层析成像重建算法;提出了绕射层析成像的非线性迭代反演法等。尽管人们对绕射层析的发展给予了极高的重视 ,但由于射线层析发展较早,方法相对成熟,因而在地球物理层析领域仍占据主导地位,而绕射层析在理论上 ,尤其在实用上还有很长的路要走。2. 国内外研究历史和发展现状井间地震层析成像技术及应用研究_郭彦刚.1917年，奠定了层析成像理论的数学基础是 J.Radon 发表了著名论文关于由函数沿某些流形的积分确定该函数（即 Radon 变换）。在此之后的几十年间，物理学、医学和天文学等领域先后对其相关领域的目标进行图像重建利用层析成像的原理，1956 年，R. N. Bracewll 在射电天文学中实现了第一个图像重建通过对图像重建二维 Fourier 反演方法问题研究，1963 年，A. M. Cormackl 发展了投影重建图像的解析数学方法利用的是 X 射线，从而为后面研制医疗 CT 设备奠定了基础，其发表了论文在发射学中函数的线性积分图像重建的应用。1971 年 G. N.Ramackandran 等研究了利用褶积替代傅氏变换的方法在 X 射线照相和电子显微照相中的应用。20 世纪 60 年代末起，在油气藏勘探开发层析技术开始得到了应用。1972 年，美国地球物理第一次对在油田上井间地震的层析实验进行了报道，这次试验并未显示出井间地震的应用潜力主要原因在于现有的实验设备简单，成像质量很差，所以没有得到业界的重视。1979 年，G. N. Hounsfield 与 A. M. Cormack 由于将层析成像原理应用于 X 射线脑颅 CT 扫描，获得了医学诺贝尔奖。层析成像得到了飞跃发展广泛的应用在各个领域。1983年，美国的许多石油公司和大学研究机构投入大量精力研究井间地震波层析成像，其动力来自于 G. McMachan 提出了井间地震波层析成像的方法和概念。相关井间地震新的算法、技术、设备等相继出现，推动力是计算机领域技术和地学领域的迅速发展，进而很多新的应用成果出现了。开始出现了许多井间地震的相关论文在勘探地球物理学家 SEG 等国际年会上。80 年代后期，斯坦福大学的一个研究小组建立了井间高分辨率反射波成像方法通过吸取地面地震反射法和 VSP 的原理和数字处理方法。八十年代，地震层析成像引起了地球物理学界学者的高度重视，自从在亚特兰大（Atlanta）召开的第 45 届地球物理勘探学家学会（SEG（Society of ExplorationGeophysicists）年会上设置了地震层析成像研究内容的专题之后，以 Daily（1984），Somerstein（1984），Pratt and Worthington（1984），Bishop（1985）等人的研究为代表，利用人工地震发射与接收系统的地震层析成像理论、方法和技术以数值模拟的形式得到深入、广泛的研究。在实际应用方面，井间地震层析成像不断扩展，在资源勘探、工程勘探、环境保护、文物调查、防灾减灾等许多应用领域都得到实验性研究并取得有效的进展，如矿井爆破（Gustavsson et al.，1986），煤层断层检测（Mason，1981；Bodoky，1985），煤矿应力监测（Kormendi et al.，1986），空洞检测（Vazqnez etal.，1985），盐丘调查（Peterson et al.，1985），坝体研究（Cottin et al.，1986）以及与核废料有关的岩体研究（Ivansson，1985；Paulsson，1985；Gustavsson，1986）等。自 90 年代开始，井间地震的得到了广泛的应用推广，基于井间地震的方法原理基本完善了，各种井中地震设备的发明为起了巨大的推动作用；包括井中震源以及使用多级检波器串的，野外采集成本大大降低了，数据采集周期得到了缩短；相应的处理软件的也相继开发出来，推动了井间地震技术迅速发展。井间地震已实际应用在油气开发与勘探领域，并得到了很好的应用效果。1985 年，Paulsson 等成功的监测了小尺度井间实验中的岩石参数，其实验的井间的距离仅有几米。1986 年，Gustavsson 在瑞典北部 Kiruna 铁矿将人工地震CT 应用于矿产探测。T.Inazzhi 等人则评价岩体的质量通过层析成像方法。1986年，McCann 等对土木工程中井间测量的应用进行了研究。针对高分辨率 CT 在1992年第62届美国勘探地球物理学会年会上开始设专题研究。当前，美国的井间地震已进入稳步发展的阶段，井间地震已从简单的二维逐步走向三维、四维(时延井间地震)，并已形成一套完整成熟的工作流程，可多井同步检测，最大井距可达 800m，可在高温井、斜井、水平井及油管井接收，并由陆上发展到海上。20 世纪 90 年代初，井间层析地震理论与技术由石油大学吴律教授首先引进国内，由他所主持的研究小组取得了许多成果在利用射线理论与波动理论进行层析反演图像重建方面。九五攻关成果“高精度弹性波与电性层析成像方法”中国科学院地质与地球物理研究所进一步发展了变速度弹性波克希霍夫积分偏移技术、基于全波信息的波动方程介质物性参数反演新方法、新技术、单分量弹性波有限元偏移技术以及井间透射波和反射波层析成像方法技术。“十五”期间，由中国石油大学、南京物探研究所和中国石化胜利油田把井间地震作为其中一个研究专题进行联合攻关研究，我国组织了“大中型油气田勘探开发关键技术研究”的国家级重点科研攻关项目。目前井间地震最主要的应用是油气田勘探开发。Lines，Harris，Tura等人利用井间地震资料的高分辨率特性来识别尖灭带,检测裂隙以及使深部的薄储集层成像;应用层析速度 井间地震资料的反射成像及槽波等研究储集层的连通性。Greaves和Paullsrn利用井间地震进行浅层蒸气驱油藏描述。近几年来，井间地震技术的独特优势己受到工程物探界的重视,国内外都有一些在地层浅部结构和地质灾害中的应用实例。如Hanming Gu等利用井间层析成像技术确定石灰岩中溶洞的形态特征。在国内,曾经有许多学者和工程技术人员作过井间地震试验。试图将井间层析成像应用到工程领域,应用于浅部地层构造和地质灾害调查。杨勤海(2005)在唐山市体育场利用井间地震波CT技术对岩溶塌陷地质灾害进行了勘察,证明井间地震波CT技术在划分岩溶,确定岩溶或岩体破碎带位置等方面是一种高分辨率的有效探测手段。王建军 曹建伟(2008)应用井间地震CT对湖北京珠国道水河大桥河中桥墩基础处基岩地质构造进行了勘察,准确地确定了该桥墩基础处断裂的位置和基本情况。国内从上世纪八十年代开始井间地震波层析成像研究，早期在反演成像算法（杨文采，1989），软件系统（朱介寿等，1994），实际应用（宋正宗等，1995）等方面取得了新的研究成果。随着研究应用的不断发展和深入，井间地震波层析方法在国内工程勘查中出现了许多成功的应用，并成功有效地应用于施工质量检测（雷宛等，2005）等方面。但井地震在隐伏断层上断点探测中的应用目前还没有相关的报道,反演方法不存在问题,其关键性难点足在松散地层的井场条件下,很难采集到主频高 重复性好的井中地震资料。总体上说，国内井间地震层析成像技术已经在油气勘探、采煤、水利工程、桥隧工程、建筑业等领域发挥了独特的作用，高精度、低成本、无损伤的优点得到了充分显示。井间地震已经逐步完成了从科学研究向实际生产力的转换。但是，井间地震层析成像技术还存在某些本身固有的缺点作为一种新型物探方法，如：有限观测角造成的层析反演多解问题以及两井之间工作距离等问题，这些问题限制了井间地震层析成像技术的发展，但也同时促进了新技术与新方法的研究。3 层析方法存在的主要问题及解决途径 成谷等, 地震层析成像中存在的主要问题及应对策略. 地球物理学进展, 2003. 18(3): 第512-518页.层析成像最早应用于医学,利用X 射线研究人体内部结构,医学上应用的巨大成功使层析成像逐渐拓展到其它领域.20世纪80年代地球物理领域引入了层析成像.地球物理领域的层析成像研究如何利用地震波、电磁波或其它场的数据对地球内部成像,但与医学层析不同,医学层析中可以获得各种角度的投影数据,而地球物理层析的发射源及接收器只能分布在地面或井中,且数目有限,因而获得的投影数据的角度范围有限. 而根据Radon 变换,只有在已知所有入射角的投影函数时才可以唯一地恢复图象函数. 正是这种先天不足(投影角度严重不全) ,限制了层析成像在应用地球物理领域的发展. 但与其它方法相比,层析成像仍不失为一种研究地下介质构造和物性分布的好方法. 层析成像在地震中的应用主要存在的问题是:(1) 数据投影角度不全; (2) 数据不均匀覆盖; (3) 反射层析中界面深度和速度的耦合; (4) 参考模型引起的非线性误差; (5) 边界作用;(6) 数据的不确定性问题(即数据中含有噪音) ; (7) 模型参数中不同分量的不确定性不同.（8） 正演方法逼近观测数据的精度问题（9）反射层析中速度与反射界面深度之间的多解性问题3 问题产生的原因和应对策略3. 1 数据投影角度不全的问题3. 1. 1 如果数据的投影角度不全,就不能唯一地恢复图象函数. 在地震层析成像中,由于震源和检波器只能离散地分布在地面或井中,且数目有限,因而获得的投影数据的角度范围有限. 决定了得到的解不唯一，只能依靠工作人员根据对工区实际情况的了解以及工作经验,施加约束来得到唯一解. 而且射线的覆盖角度有限使得层析反演出的图象沿着射线路径模糊. 地震勘探的观测方式决定了数据投影角度不全的程度,对于地震勘探中的三种典型观测方式VSP(垂直地震剖面) 、Crosshole (井间地震勘探) 和SRP(地表反射地震勘探) :Crosshole 投影角度一般在50 度之内;VSP 投影角居中;而SRP 投影角最小,仅为25 度.地震勘探各种观测方式中数据的投影方向决定了地震速度层析成像在各方向的分辨率. 在反射地震勘探中,由于射线主要是垂直方向的,因此在反射地震速度层析成像中,垂向分辨率不如横向分辨率;而井间地震勘探中,由于射线主要是水平方向的,因而横向分辨率不如垂向分辨率. 即反射地震勘探存在垂向模糊性(Vertical Smearing) ,而井间地震勘探中存在横向模糊性(Horizontal Smearing) 1 . 反射地震勘探中的射线多为垂直方向的,因而能较好地确定了地质体垂向部分的左右边界,即地质体垂向部分左右边界清晰,而由于水平方向射线较少,因而地质体水平部分的上下边界相对模糊,即垂向模糊. 由于多数地质体都是水平方向分布的,因而在反射地震勘探中层析反演出来的图象将深受这种“垂向模糊性”的影响. 同样的道理,由于井间地震勘探中水平方向的射线居多,因而地质体的水平部分的上下边界能够较精确地确定,而垂直部分的左右边界相对模糊,即横向模糊.3. 1. 2 在实际应用中,主要从采集和处理两方面采取措施来减少数据投影角度不全的影响. 在采集方面可以通过增大检波器排列的长度增大投影角度;在数据处理方面发展更多扩展空间谱范围的方法,如补零法、对称扩展法、G2P 外推法和最大熵逼近法等,这些方法在某种程度上能够增大数据的投影角度,但作用有限. 此外多种观测方式获得的数据联合起来进行层析也可以增加数据的投影角度.3. 2 数据不均匀覆盖问题3. 2. 1 不均匀的射线覆盖主要有两方面原因:一方面与观测方式和数据采集的质量有关,如炮点之间不等距、检波点之间不等距、不同炮的数据偏移距不同、不同道的数据之间质量不同等. 另一方面,对每一炮而言,所有的射线都从炮点出发,因此在炮点处射线是非常密集的,而在远离炮点的地方,射线相对稀疏. 此外速度结构也影响射线的覆盖方式,由于速度的非均匀性射线可能在某些地方聚集而在其它某些地方发散2 . 在地震勘探的观测方式中,SRP 中射线的覆盖是相对均匀的,而VSP 中射线覆盖最不均匀,因此VSP 层析受数据不均匀覆盖的影响最为严重.在层析反演中,如果对数据的不均匀覆盖不进行处理,则不均匀覆盖的数据会以一种不均匀且不可预测的方式使反演结果发生畸变2 . 数据的不均匀覆盖表明数据对不同模型参数的控制程度不同.而层析成像中最常用的常规DLS(阻尼最小二乘) 法的最优准则是使数据拟合差和解长度同时最小. 而常规的DLS 法中解长度的定义为mm ,即所有的模型参数具有相同的权重,因而常规的DLS 法忽略了数据的不均匀覆盖的影响,Wang and Braile2 证明用忽略数据不均匀覆盖影响的DLS 法得出得模型校正量正比于射线的覆盖程度. 射线覆盖程度高的区域内的网格点被过分校正,而射线覆盖程度低的区域内的网格点校正量不足,而那些没有任何射线穿过的网格点的校正量为零.3. 2. 2 减少数据不均匀覆盖对层析结果的不利影响的途径是减少模型校正量和射线覆盖程度的相关性. 主要有两种方法,一种是根据数据的控制程度对模型参数进行加权,如变阻尼法2 . 另一种是采取某种措施使数据对每一模型参数的控制程度近似相同,如变网格采样法2 . 变阻尼是应用不均匀的先验信息来匹配不均匀的数据分布,而变网格采样法是重新调整数据信息的分布使数据对每一参数模型的控制程度更平均. 二者的选择标准是基于对射线覆盖程度的定量分析. 对于某些网格点没有任何射线穿过的情况,变阻尼法本身具有不稳定性问题. 有两种方式可以来增加变阻尼法的稳定性:设定一个最小阻尼;平滑法和变阻尼法联合应用.3. 3 反射层析中界面深度和界面上方速度的耦合问题3. 3. 1在反射地震勘探中,纪录的数据受两个过程的影响3 :通过界面上方速度场的透射过程和在界面处的反射过程. 这两个过程对数据的影响耦合在一起,很不容易分离. 典型的例子就是我们不能确定零偏剖面上的旅行时变化到底是由界面上方的速度变化引起的,还是反射界面深度的变化引起的.速度和其下的反射界面的耦合程度有强有弱.当界面上方只有小尺度的速度变化时,速度和反射界面的耦合程度较弱,而当界面上方有大尺度的速度变化时,二者的耦合程度较强. 此外速度变化与反界面之间的距离也决定了速度和反射界面之间耦合程度的强弱. 一般来说,速度变化越接近反射界面,二者之间的耦合程度越强. 而层析成像的目的是为偏移提供精确的速度输入并将界面深度和界面上方的速度解耦. 因此有必要研究速度和反射界面解耦的分辨率.理论上,反射面位置引起的旅行时变化正比于入射角(入射射线与垂向所成的角度) 的余弦值,而速度变化引起的旅行时变化正比于入射角余弦的倒数,速度变化和反射面位置是可以解耦的. 但因为速度和反射面位置引起的旅行时变化随入射角(偏移距) 变化很微妙,而且由于其它因素的影响(如数据覆盖角度不全、边界作用、非线性误差的影响等) ,实际上速度和反射面很难完全解耦. 其解耦的分辨率取决于构造形态、地质约束、离工区边界的距离、信号的质量、数据的密度,以及数据的最大偏移距4 .其中构造形态决定模型中速度变化的尺度、速度异常的形状和幅度、与反射界面之间的距离等(一般速度和反射界面深度的解耦能力随速度变化在深度上的增加而降低) ;地质约束决定补偿数据投影不全问题引起的解非唯一性的先验约束;离工区边界的距离决定受边界作用影响的程度;信号的质量决定受数据中包含的噪音的类型和大小的影响;数据密度决定数据对参数的控制程度;最大偏移距决定数据的最大投影角度.3. 3. 2本征分析是研究速度和反射面解耦的分辨率的有效工具. 本征分析4 6 将速度和反射面之间相互影响的关系分为:速度和反射面之间建设性地干涉作用(两者之间的影响是积极的,有助于恢复正确的速度和反射面位置,但只限于恢复趋势;对应较大的本征值,在相对少的迭代次数内即可反演得到) 、速度和反射面之间破坏性的干涉作用(两者之间的影响是消极的,对应非常小的本征值,迭代反演中很难收敛,在反演结果中产生假象) 和速度变化对反射面只有很小的影响(两者之间是不太相干的,彼此互相影响的程度不大,本征值大小介于上面两种情况之间) .界面深度和界面上的速度分别对应地下介质的高、低波数信息,二者互相补充但不能互相代替. 大多数的层析方法一般都能够比较好的解决其中一个问题,利用透射能量可以求解大尺度的速度场,利用散射层析(在地面反射地震勘探中,偏移就是散射层析的一种形式) 可以求解小尺度的速度场. 这里的大尺度、小尺度以地震数据的菲涅尔带范围作为分界.但这里隐含了一个基本假设,即如果致力于求解一种信息,则另一种信息必须是相对精确已知的. 界面深度和速度二者之间的互补性使得在分立求解方法中二者互为输入,且无论求解哪一种信息的方法均对输入的微小误差十分敏感. 具体地说,要求准反射界面的位置,精确的速度模型是必需的,而要求得精确的速度场,又要求反射界面的位置十分准确. 地震数据中透射作用和反射作用内在的耦合使应用层析方法进行叠前速度分析十分复杂化.针对界面深度和速度二者之间的互补性,解决这种耦合作用的方法也主要分为两种:由于求解高波数信息的方法和求解低波数信息的方法均对输入中的误差非常敏感,因此当界面上方只有小尺度的速度变化(意味着当进行高波数信息反演时,低波数信息的误差不算太大) 时,可以采取低波数信息(速度场) 和高波数信息(界面深度)迭代求解,由透射层析法求解速度场(界面位置作为输入) ,散射层析法(主要指偏移) 求解界面位置(速度场作为输入) . 而Stork 和Clayton (1985) 7 也证明了反射界面小尺度的变化虽然会在层析反演中引起假象,但这些假象或者很容易识别,或者对偏移只产生小的影响. 表明层析方法和偏移迭代应用解决小尺度速度变化和反射界面弱耦合的情况. 尽管如此,由于求解速度场和求解界面位置的方法均对输入的误差非常敏感,这种先固定一个,再求取另一个的迭代解法很难在合理的次数内收敛.而当界面上方有大尺度的速度变化时(意味着当进行高波数信息反演时,低波数信息的误差不能算小) ,这时速度和界面深度耦合程度很强8 ,因此只能采取界面位置和速度同时反演的方法. 这时是一种具有速度和反射面深度两类不同参数的非线性反问题. 研究证明当非线性反问题涉及到很多不同类型的参数时,解迭代类方法最常用的梯度方法把具有不同特征和物理量纲的参数混在一起,导致收敛性很差,并且算法的收敛性对不同参数类型之间的规一化( scaling) 具有很强的依赖性. 因此本文推荐采取子空间法9 . 子空间是由在模型空间内的基矢量7 ,8 张成的空间,该方法可以灵活处理具有相同类型参数和不同参数类型的大型反问题,均衡考虑目标函数对不同参数类型的不同的依赖性,当前模型的扰动量并不依赖于单个模型参数类的规一化(scaling) . 子空间方法是灵活的,可视为梯度法和矩阵法的一种交叉. 子空间方法也是一种很好的参数化方式,因为该方法通过把参数表示为基矢量的组合形式,将求大量的参数的问题转化为求有限的基矢量的组合系数的问题,大大降低问题的维度,使问_题的维度与参数类型的个数同级. 当然,在子空间方法中,基矢量的选择是至关重要的. 基矢量为目标函数对每一类模型参数的梯度,如果模型参数类小于4 ,则可以加入梯度变化率的信息作为额外的基矢量,使得对于每一种参数类型至少有一个基矢量与之有关,其具体构成具体见原文. 除了在参数化方面和处理不同类型参数方面的优点外,该方法允许在不同参数类型之间以及同一类型参数的不同分量之间具有相关性,表现在模型协方差矩阵的构成上. 这在同时反演P 波速度、S 波速度和密度的方法中非常有用(因为P 波和S 波速度是相关的,在先验模型中包括这一点将对算法的收敛性和层析结果的精度大有裨益) .但这种方法也具有所有的全局寻优方法的缺点,即不能保证得到的解为全局最小解.3. 4 参考模型引起的非线性误差问题3. 4. 1 实际应用中由于不知道正确的射线路径,解决的最佳方法是在最佳的猜测模型(参考模型) 上进行射线追踪,因此参考射线路径和真实射线路径之间存在差别,这一差别在反演结果中引入了误差. 参考模型引起的误差是一种非线性误差,它影响反演所用的线性系统. 参考模型引起的非线性误差不同于线性系统的不完整性引起的线性误差,这二者是可以区分的.在反射层析中,参考模型中包括速度和反射面位置两类参数. 而参考模型中速度和反射面位置的不正确引起的非线性误差对反演结果的影响是不同的4 . 参考速度与真实速度的偏差对射线路径只产生较小的影响,只会使反演的图象模糊;而参考反射面与真实反射面位置的偏差会严重影响射线路径,对反演结果的影响不可预料. 参考速度引起的非线性误差的大小与速度变化偏离震源和检波点的位置有关. 当速度变化离震源和检波点的位置较远时,射线更可能改变它的路径. 因此速度变化偏离震源和检波器的位置越远,真射线和参考射线之间的偏差越大,非线性作用越强. 而参考模型中反射面位置引起的非线性误差大小与参考反射面的平滑与否有关,如果参考模型中的反射面位置是平滑的,则引起的非线性误差相对较小,而如果反射面位置是不平滑的,则引起的非线性误差较强4 .地面反射地震和井中地震两种观测方式相比较,前者受参考模型引起的非线性误差的影响最小.而井间地震层析中最严重的问题是不正确的射线路径引起的非线性误差问题. 因为地质体大多是水平走向的,而在地面反射地震勘探中射线多为垂直方向的,故地质体的垂直方向上的变化对射线路径误差的影响较小,而地质体水平方向上的变化对射线路径误差的影响较大. 而在井间地震勘探中,由于射线多为水平方向的,情形正好相反.3. 4. 2 改进射线路径误差引起的非线性误差的方法:表现在模型的选取上是利用尽可能多的先验信息,选择最接近真实模型的参考模型;而表现在算法上是在改进后的速度模型上重新进行射线追踪和反演,即用迭代线性反演来逼近真实的非线性反演.3. 5 边界作用边界作用是由于在边界处数据的覆盖程度不足引起的,并由于大偏移距射线的存在而使边界作用影响工区的中间区域,使远离边界很大一段距离内的图象都发生畸变,对反演结果产生非常不利的影响. 减少边界作用的办法:一是增加探区宽度(增大检波器的排列长度) ;二是在探区的边界施加一些约束.3. 6 数据的不确定性问题3. 6. 1 测量的数据不可能是很精确的,测量误差的存在或其它因素的影响,使我们获得的数据中含有不同类型或不同能量的噪音. 这些噪音在反演过程中会产生很大的影响,尤其是在小本征值反演中. 含噪音数据的层析反演会使层析结果中包含很多起伏,这些起伏没有地质意义,除了淹没真实的构造形态外,这些高波数的起伏也会产生屏蔽区和焦散区,使以后的迭代过程(涉及到射线追踪) 不稳定. 在层析反演过程中,处理数据中的噪音(或称数据的不确定性) 一般是以数据协方差矩阵的形式根据数据的方差对数据进行加权. 但这个数据协方差矩阵很难精确确定,因此在层析反演过程中的应用有一定难度.3. 6. 2笔者推荐两种方法:一种方法表现在未知参数的处理上10 ,11 ,把数据中的噪音当作一种静校正量作为另一种未知参数,同时反演速度和静校正量.这是一种很有价值的处理方式,该方法的优点是:由于把静校正项(噪音) 作为未知参数矢量,因而认为数据中是没有误差的,避开了确定数据协方差矩阵的难度,在层析反演速度的同时,又可以得出静校正量的估计值; 另一种方法是对目标函数的优化上12 . 考虑地震数据中包含的噪音类型可能不同,因此在目标函数的优化上作了改进,使目标函数联合l1/ l2 范数最小,该种方法基于l2 范数适合于处理高斯型分布的噪音,而l1 范数适合于处理具有大误差噪音数据,能针对数据中的噪音是高斯型分布的,还是存在有大误差噪音的不同情况灵活处理.3. 7 模型参数不同分量的不确定性不同的问题3. 7. 1 在层析问题中,尽管我们得到的系统是超定系统,但对于模型的解估计来说却可能是一个混定问题,即对于模型中有些分量是超定的,有些分量是欠定的,有些分量是处于零空间的. 也就是说,模型的不同分量其不确定性是不同的. 超定的分量不确定性弱,欠定的分量不确定性强,而对于处于零空间的模型分量不确定性非常强. 前面所提到的各种因素都对层析反演中模型参数的不确定性问题有贡献. 对前面各种因素处理得好,则模型参数处于零空间和欠定的分量就少,而如果处理得不好,则模型参数处于零空间和欠定的分量就多.一般来说,模型中不同分量的不确定性要有一个先验地估计. 要使模型的解估计与真实模型不致偏差太大,必须对模型参数中欠定的及处于零空间的分量进行约束. 具体问题的欠定谱分量和零空间谱分量通过本征分析来大致确定. 在本征分析中,对应大本征值的模型参数矢量基本上是超定的,在迭代反演过程中很快收敛;对应小本征值的模型参数矢量基本上是欠定的分量,在迭代反演过程中收敛很慢;而对应零本征值的模型参数矢量对应着处于模型零空间的分量,这些分量对数据没有影响,因而是无法反演出来的. 在很多最小二乘反问题中,零空间分量和欠定分量集中分布在某个谱带之内13 ,如在井间层析中零空间和欠定谱分量主要与水平波数有关,而在上面提到的速度和静校正量同时反演的问题中,当用直射线和井间观测方式时,对慢度扰动参数来说,水平波数是欠定的,而对静校正量参数来说,直流分量位于零空间,对于某些其它模型来说,静校正量的低波数分量也是欠定的10 . 在波形反演中,甚低波数的速度变化处于零空间.模型参数中欠定的分量和处于零空间的分量有时对于解估计逼近真实模型很重要,因此在层析反演过程中,对这些模型参数中欠定的和处于零空间的分量通过约束补充以减少这些分量的不确定性十分必要. 不进行特别处理,轻则在层析反演结果中引入假象,重则使层析反演结果完全畸变.3. 7. 2 解决模型参数不同分量之间的不确定性不同的问题的方法是加入平滑和规则化约束.平滑方法主要分为静态平滑法和动态平滑法.静态平滑法指在层析反演的过程中平滑因子保持不变,而动态平滑则在层析反演的过程中根据具体况调整平滑因子的大小和形状. 动态平滑14 ,15 的效_果要优于静态平滑. 动态平滑方法有很多种,如多重网格法14 和褶积压制法15 ,多重网格法是一种收敛速度很快且计算量较小的方法,其收敛速度是线性的且与网格大小无关,而很多其它的迭代方法的收敛速度均依赖于网格大小且随着网格大小趋近于零而恶化. 其计算量小是因为该方法只在最粗的网格上进行精确计算,而在细网格上都通过下落到粗网格上近似计算. 该方法首先从细网格下落到粗网格,然后又重新上升到粗网格,从细网格到粗网格由限制映射实现,由粗网格到细网格由延拓映射实现.整个过程基本上由前光滑、粗网格校正和后光滑三个步骤实现16 . 而在褶积压制法中平滑子作用(一种褶积作用) 于射线路径矩阵的每一行,使射线路径被拓宽成一个带,带的宽度由平滑函数的宽度决定,即褶积压制法认为邻近网格点的速度值是高度相关的,该方法实际上是控制分辨率和方差之间的折衷程度的一种方法.规则化方法一般分为阻尼型规则化和导数型规则化. 当在层析反演中采取阻尼性规则化方式时,普通的CG(共轭梯度) 法很容易收敛,而当在层析反演中采用导数型规则化方式时,CG 法不容易收敛.VanDecar 和Snieder (1994) 对采用导数型规则化方式的CG法进行了改进,对CG法的搜索方向进行预处理,并利用规则化方程中的先验“构造”信息,使该算法只依赖于模型参数的个数而与数据的个数无关,大大减少了迭代次数17 . 采用导数型规则化方式得到的解为最少“构造”解,与在模型参数的先验估计不能比较精确地给出时按照最接近初始模型的选择标准产生的解估计(这样的解估计中包含的“构造”量比解释数据必需的“构造”量更多,称为更多“构造”解) 相比较,最少“构造”解中只包含能够与数据匹配的最少的“构造”量(最不复杂的解) 17 . 模型中包含的“构造”量由一阶或二阶导数定义. 一阶导数约束得到最平解,二阶导数约束得到最平滑解.Phillips 和Fehler (1991) 针对线性问题对很多常用的平滑和规则化方法(中值平滑法、空间平均平滑法、褶积压制法、网格移动平滑法、一阶差商规则化法、仅在水平方向一阶差商规则化法、IRIS 法、DLS法等使用单约束参数的方法) 的效果进行了比较18 ,认为一阶差商规则化法在逼近真解的程度和处理不同噪音数据的稳定性两方面表现最好,且对高噪音数据和低噪音数据都适用. 褶积压制法是以DLS 法为基础的方法中在逼近真解的程度和稳定性两方面表现最好的,但平滑的形状限制了该方法能_实现的分辨率的上限. 而空间平均平滑和中值平滑各有特点,空间平均平滑后的图象较平滑(“构造”更少) ,而中值平滑后的图象“构造”更多些.在实际应用中,模型参数的不确定性一般由模型协方差矩阵表示,即各种平滑和规则化约束主要体现在模型协方差矩阵的构造中. 同数据协方差矩阵类似,模型协方差矩阵也很难精确确定. 应用中使用的模型协方差矩阵大多数是对角的,这意味着模型的各个参数之间是不相关的,但实际情况并非总是如此,很多情况下不同类型参数之间是相关的,同一类型的参数之间也可能相关. 因此使用非对角的模型协方差矩阵十分必要,但这会在实际计算上造成困难,因为实际计算要求模型协方差矩阵的逆. 本文推荐一种构成非对角模型协方差矩阵的十分实用的方法,模型协方差矩阵的逆可以作为线性滤波器矩阵的自协方差得出,而避免了直接计算模型协方差矩阵的逆19 ,滤波器对应着不稳定的谱分量,即只对处于零空间或欠定的谱分量进行约束. 通过把先验信息的应用解释为一个简单的线性滤波运算,则构造非对角的模型协方差矩阵的问题归结为构造零空间和欠定分量的傅立叶系数矩阵的问题. 一旦理解了这种线性滤波器的解释,则对于任何的病态反问题构造非对角阻尼矩阵都是很直接的.3.8 正演方法逼近观测数据的精度问题在层析问题中,观测数据与正演模拟数据的拟合程度是评价解估计的一个重要准则,因此正演方法逼近波在地下传播的真实过程的程度就成为关键。逼近波真实过程的程度越高 ,则求得的解估计中由正演造成的误差就越小。在绕射层析中 ,应使用尽可能符合真实地下介质情况的波动方程。这也反映了绕射层析的一个重要的发展方向: 由标量波、声波绕射层析向弹性波绕射层析发展。射线层析中采用的正演方法多为射线追踪 ,因此射线追踪的精度对射线层析的影响很大 。 射线追踪的方法很多 ,要针对所采用的反演方法选择一合理的射线追踪方法 。 此外 ,射线层析中由于理论基础的限制也带来了一定的误差 。一般来说 ,地震射线层析中 , 旅行时是慢度沿着在二维或者三维空间中的一条射线上的线积分 。这简化了数学处理 ,但它与我们的物理直观知识是相矛盾的 ,事实上无限窄射线仅存在于无限小波长( 0) 的情形 。由于地震数据的带限特征 , 在地表接收到的旅行时信息不反映沿着单一射线路径的速度变化 , 而是反映在单一薄射线 ( thin ray) 周围一有限体积范围内的速度特征 。即地表接收的旅行时是沿着二维或三维空间内薄射线周围某一有限体积范围内的体积分 ,为此人们发展了波路与菲涅尔体积的概念 ,给出了计算波路和菲涅尔体积的方法 ,并将其应用于旅行时层析 , 得到了比利用薄射线层析更为满意的结果 。3.9 反射层析中速度与反射界面深度之间的多解性问题到目前为止 ,针对反射层析中反射界面的处理主要有 3 种方式 : 假设在层析反演速度的过程中反射界面的深度是精确已知的 ,层析反演只是迭代求解速度的修改量 ; 把反射界面的深度也作为未知量和速度一起一步反演 ; 通过偏移和层析迭代求解速度和反射界面深度 ,在求取速度的修正量前假设反射界面的深度已知 ,在用层析方法求出更为精确的层速度之后 ,用新求出的层速度再利用偏移求取更精确的反射界面的深度 ,如此迭代进行直到满意为止 。其中第 1 种处理方式由于误差较大 ,近来已很少使用 ; 第 2 种处理方式在解决速度和反射界面深度的多解性方面较第 3 种方式差 ; 第 3 种处理方式因为涉及到层析求取速度修正量和偏移求取反射界面深度的交互进行 ,尽管与前两种处理方式相比更精确 ,但计算量较大 。3.10讨论和结论地震层析成像中存在的这些问题各有特点,但又互相联系. 它们可以使层析反问题的性状变坏. 由于反问题性状的好坏是影响反问题稳定性的重要原因,因此在层析反演的前期及反演过程中针对这些因素进行处理非常必要.3. 10.1 增大检波器排列的长度、发展多种扩展空间谱范围的方法及多种观测方式获得的数据联合层析.数据投影角度不全问题.3.10. 2 根据数据的控制程度对模型参数进行加权和采取措施使数据对每一模型参数的控制程度近似相同的方法,减少数据不均匀覆盖对层析结果的不利影响.3.10. 3 反射层析中速度和反射面位置之间的耦合性是由于地面反射地震勘探中透射作用和反射作用的耦合产生的. 针对这种耦合性主要有两种处理方式:一是速度和界面深度迭代求取的方法,适合于界面上方只有小尺度的速度变化的情况,但所需迭代次