2001年9月地球物理学进展第16卷第3期井间地震层析成像的现状与进展裴正林（石油大学（北京）物探重点实验室，北京，100083）摘要：综述了井间地震层析成像研究的现状，给出了小波变换域井间地震层析成像方法的最新进展，并对井间地震层析成像研究给予展望.关键词：井间地震层析成像；小波多尺度；研究进展中图分类号：P315.3+1文献标识码：A文章编号：1004-2903（2001）03-0091-071井间地震层析成像的研究现状井间地震层析成像也称为井间地震CT技术，它能够提供被探测地质体的构造和岩性分布的高分率图像.井间地震CT技术是从医学CT技术发展起来的，其数学基础是Radon变换.井间地震CT的研究基本始于20世纪70年代初，80年代处于对大量模型数据和少量实际数据的成像研究阶段，90年代以来，井间地震CT进入实用化阶段，并取得不少可喜成果，同时，也逐渐意识到射线CT所固有的缺点，开始研究波动方程CT.从地震波的运动学和动力学特征出发，井间地震CT方法可分为两大类：一类是基于几何光学或射线方程的方法称之为射线CT；另一类是基于波动方程的方法称之为波形CT.当非均匀体的线性尺度大于地震波长时，射线CT是适用的；而当非均匀体的线性尺度与波长相近时，衍射和散射就起主导作用了，基于射线理论的成像方法就不再适用，这时候必须用波动方程CT方法.井间地震层析成像方法主要包括两部分：正演方法和反演方法.井间地震层析成像的正演方法可分为两种；一是射线追踪方法；二是波场的数值模拟方法.射线理论和射线方法是研究地震波传播理论的重要方面之一.用射线理论可以研究地下复杂构造、横向不均匀介质中的地震波传播问题.经过射线追踪，计算地震波的走时、波前和射线路径.70年代以前的各种射线追踪方法一般适合于较为简单模型的射线追踪［1］.由于实际的介质速度变化较大（速度差大于10%），因此，需要研究复杂结构模型的射线追踪方法.收稿日期：2001-03-15；修订日期：2001-06-15.基金来源：“九五”国家科技攻关项目资助（959130602）.作者简介：裴正林，1962年生，2000年获中国地质大学（北京）地球探测与信息技术专业博士.高级工程师，现在石油大学（北京）从事博士后研究.主要研究方向：信号处理，小波变换、遗传算法及神经网络应用，层析成像理论方法和地震数据处理、偏移方法等方面研究.E-mail：zhenglinpei@.近年来，这方面的研究取得新的进展.目前应用效果好的弯曲射线追踪方法有两类：一类是基于程函方程的有限差分方法；另一类是基于Huygens 原理的最短路径射线追踪法.基于程函方程的有限差分方法有有限差分方法的E-O 格式［2，3］，它能以较高的精度计算较复杂模型的初至波走时.线性插值射线追踪法（LTI ）［4，5，6］是程函方程有限差分射线追踪［7，8］的改进.这类方法是全局算法，能模拟较复杂介质，编程简单、计算速度快、精度较高.基于Huggens 原理的射线追踪有两种全局算法：一是基于网络理论的最短路径算法［9，10］；二是基于动力学的波阵面算法［11，12］.这两种算法都能模拟复杂介质，而且是一次计算即可得到全部网格节点上的最短路径和最小走时，计算效率高.井间地震波场的正演数值模拟是求解变系数偏微分方程的过程.目前，解偏微分方程的数值方法中较常用的有：有限差分法，有限元法，边界元法，伪谱法，积分法等［13］.其中有限差分法比较容易程序实现，并且具有较快计算速度和较高计算精度，很适合在波动方程层析成像中使用.目前，井间地震层析成像2-D 正演方法基本趋于成熟.井间地震层析成像是非线性反问题.井间地震射线CT 是通过逐次线性化反演来实现的，即可归结为求解一个大型的、稀疏的、常常是病态的线性方程组.为了在这种情况下求出满意的解来，必须从数学上把病态问题转化为良态问题求解，并把更多的先验信息引入到反演过程中.近十年来，许多研究者不断研究发展新的方法，改善了在数据不完全，有噪声和初始模型不准情况下的成像效果.射线CT 从直射线层析成像方法［14］发展到弯曲射线层析成像方法［10，15，16］；反演万法由最小二乘法发展到各种约束条件下的阻尼最小二乘法［17，10，15，18］以及统计法如最大熵法［19］；观测参数数据也由单一的走时向多种参数数据（P 波、S 波走时、振幅、信号谱等）发展.井间地震P 波、S 波走时层析成像方法［20］以及多震相成像方法如直射线的反射十透射井间层析成像方法［21］及弯曲射线的反射十透射井间层析成像（CCRTT ）方法［22］；井间地震和逆VSP 联合层析成像方法，一定程度上克服了井间地震透射层析成像的多解性及有效成像区域小的缺陷［23］.衍射CT 是波动方程CT 的一次近似（假定弱散射，ΔV <15%的小扰动）.有两种方法：Born 法和Rytov 法.Born 近似是利用波场的振幅，而Rytov 近似则利用波场的相位.与射线CT 相比，衍射CT 提高了分辨率，可以减少由于有限观测角所造成的假像［24，25，26］.在一定条件约束下，Born 近似和Rytov 近似的效果是相同的；在高频假定了，Rytov 近似层析与走时层析是等价的［27］.声波衍射CT 方法可以获得纵波速度的像，或密度和压缩系数的像［28］.弹性波衍射成像方法可以获得纵横渡速度或密度和拉梅常数的像［29，30，31］.目前衍射CT 大多只用于均匀的背景介质，而变背景的比较少［32，33，34，35］.散射CT 方法是广义Radon 变换的一级近似逆.它可用于波还不均匀的背景介质.散射CT 方法在己知准确的背景速度分布时可以给出分辨率很高的成像结果，但在不能给出准确背景速度场时效果很差［36，37］.为了得到准确的背景速度分布，可以采用级联算法［38］.实际应用时，衍射CT 和散射CT 方法都存在一些有待解决的问题，如散射波的提取、对波形产生影响的各种干扰因素的消除（震源函数、介质的吸收、震源辐射、接收器的耦合）［39］，而且这些成像方法只有当散射体与背景的速度差较小时才有效，因而使其实际应用受到限制.而波形CT 方法是利用全波场信息，能够正确地反演地下的真实介质模型，因此，具有・29・地球物理学进展2001年潜在的科学价值.近年来，波形CT 方法成为研究的热点，并有所发展.波形CT 又可分声波方程CT 和弹性波（全波）方程CT.声波方程是弹性波方程的近似.对2—D 声波方程速度反演研究比较多［40，41，42，43］.虽然弹性波方程CT 给出的图像略好于声波方程的，且给出V p ，V s图像，但声波方程的计算时间约为弹性波方程计算时间的1/8［40］.由于波形反演的目标函数中存在大量的局部极小，因此，波动CT 存在着收敛速度慢，对初始模型依赖性强以及易于陷入局部极小的缺陷.有一种波动方程走时反演方法［44，45］.该方法是利用地震波的走时差来代替波形反演中的波形拟合差进行波形反演的，提高了抗噪能力和速度成像的分辨率（相对于射线走时成像而言）.另外，加快了收敛速度.但这类方法仅利用了波形的走时信息，从某种意义上说，它减少了波形反演的局部极小个数，但这是以降低波形反演的分辨率为代价.射线CT 获得是速度分布的光滑估计即速度模型的低频成分，提供低分辨率图像，而波动方程CT 是估计速度模型的不均匀性和不连续性即高频成分.因此，二者联合使用则很适合于高分辨率的图像重建［40］.井间地震层析成像问题的非线性是与介质的非均匀性联系在一起的.对于小扰动非均匀介质，常规的线性化反演方法可以得到较好的解.而对于大扰动的非均匀介质的成像问题，采用常规的线性化反演方法，反演结果不好，这是由于在此种情况下，线性化近似已经不再有效.可以采用逐次线性化反演方法，其结果尚可满意［46］.井间地震层析成像是非线性反问题，具有非线性性、多解性与不适定性、高维数（为达到分辨率要求，模型网格点数量很大）；另外，还存在弯曲射线、大扰动、有限频带以及有限观测角（观测数据不完备）等问题，因此，井间地震层析成像是复杂非线性反问题.井间地震层析成像反演方法可分为两类：基于算子的线性化反演方法和基于模型的完全非线性反演方法.线性化或拟线性化的反演方法又称为确定性反演方法如LSQR 法、梯度法、牛顿法等；完全非线性反演方法又称为随机反演方法如Monte Carlo 法、模拟退火算法（SA ）等.确定性反演方法一般是局部寻优的，而随机性反演方法多为全局寻优.线性化CT 反演方法如正则化方法，Born 和Rytov 近似，散射法等，通过忽略高波数成分即局部极小，建立了适定的算法，换言之，线性化反演方法是通过忽略参数的部分非线性成分，改善了反问题的适定性，反演出有意义的解.因此，线性化CT 反演方法只有当初始模型在目标函数的全局极小的邻域内（在真实模型附近）时，近似线性化的关系才成立，即只有当引人很好的背景速度场先验信息约束时，线性化CT 反演结果才能接近全局最优解，否则易陷入局部极小.因此，现行的井间地震层析成像方法缺点是，对初始模型的依赖性强，即图像的质量受初始模型的控制，仅适合于低对比度（速度差小）异常体成像，图像分辨率偏低以及不可缺少的先验知识等.很显然，对于复杂模型和大扰动非均匀体成像问题，现行的井间地震层析成像方法无法得到好的图像.多重网格法是目前多极值目标函数反演的有效方法之一.可以有效地改善常规一次性反演方法的性能.已有不少研究者进行了这方面的研究.地震反射波多重网格法2-D 反演［47，48，49］，井间地震走时波形联合反演方法是解决波形反演中局部极小问题另一种方法［50］.尽管多尺度方法一定程度上提高了线性化层析成像图像的分辨率和图像质量，但仍未・39・第16卷第3期裴正林：井间地震层析成像的现状与进展完全克服线性化反演方法的不足，自然无法彻底地解决非线性问题.因此，地球物理学者将随机概念引入非线性问题反演，而且取得了不小的突破，出现了许多现实的反演方法，如Monte Carlo 法、模拟退火法（SA ）、遗传算法（GA ）等.随机反演方法对目标函数不作任何数学要求（如连续性、可导性），适合于复杂的多极值目标函数的全局优化.将SA 法和GA 法应用于1-D 地震波形反演［51，52，53，54］.应用SA 法进行地震反射成像研究［55］.应用GA 法研究了2-D 地壳速度结构反射走时反演［56］.采用SA 法进行简单模型走时层析成像研究［57，58］.显见，二维井间地震层析成像的随机反演算法处于研究阶段，其应用实例也很少.对于复杂非均匀介质，现行的井间地震层析成像方法实际应用时往往图象分辨率不高、效果也不理想.2井间地震层析成像的进展2.1小波多尺度地震层析成像小波分析是当前研究非线性科学的主要工具之一.小波变换也是一种多分辨变换，它能够将信号的低频成分映射到由大尺度小波基张成的粗采样子空间，把高频成分映射到由小尺度小波基张成的较细采样子空间，因此，小波变换为多重网格法提供了一种自然框架.小波变换的这一特性已经应用于1-D 地震波形反演［59，60］.对于井间地震层析成像非线性反问题，一种方法是通过小波分解压缩速度模型，来提高成像计算速度［61］；另外一种方法是在小波域对线性化方程进行多尺度分解的小波多尺度地震层析成像方法［62］.该方法将地震层析成像反问题分解为不同尺度上的一系列反问题.大尺度的反问题的目标函数局部极小少而且互相分得很开，因而用常规的反演方法可以较容易地收敛到全局极小的邻域内；将大尺度的解当作次一级尺度反问题的初始模型，再进行反演.如此类推，直到求出对应于最小尺度的原反问题的最优解.小波变换多尺度方法有效地改善了常规线性化反演方法依赖初始数模型、易陷入局部极小的不足，是复杂介质和大扰动非均匀介质成像的有效方法，另外，小波多尺度方法还能够提供不同分辨率的反演图像，有益于图像的解释.2.2井间地震多尺度遗传进化层析成像针对井间地震层析成像反问题所具有的特点：模型参数维数高，且参数在空间分布上具有一定的规律性，模型可行解空间巨大等，把简单遗传进化算法（Evolutionary Programming ，简称EP ）［63］应用于井间地震层析成像反演时应进行改进［62］.（1）二维实数编码地震层析成像问题类似于二维图像的恢复问题.一个个体就是一幅图像.此时，若采用一维染色体编码则很不方便，尤其是变异操作很不直观.因此，采用二维实数编码.每个二维染色体就对应于一个二维速度模型，每个实数编码的基因就表示对应的速度网格节点上的速度.（2）模型扰动方法在遗传进化算法中，采用模型局部扰动方式作为对父代染色体的变异.二维染色体的变异不同于一维染色体的点位变异，而是采用基因块变异或模型速度块扰动.模型速度块大小的选取的方法有两种：随机大小速度块和固定大小速度块.随机速度块，其大小可以从一个像元到整个模型变化.模型的速度块的速度扰动范围是自适应可变的，最大范围从模型・49・地球物理学进展2001年速度的下限到上限.这样一来，染色体（速度模型）的变异自动实现约束，并约束在速度参数值的选择范围内，以提高计算可靠性及计算速度.（3）多尺度逼近遗传进化算法我们将小波多尺度分析思想也引入到随机最优化过程中，将多尺度反演和进化算法反演的结合起来，称为多尺度逼近遗传进化算法（Multiscale Evolutionary Programming Algorithm ，简称MEP ）.该方法的有三个特点：一是对原反问题分解为尺度由小到大的反问题序列.随着尺度的增大，对应反问题的解空间范围缩小，这样以来，加速了反演搜索速度，减少了该尺度反问题的求解计算成本.二是对于大尺度反问题，可以通过小尺度来细化.三是与常规遗传进化算法所采用的固定搜索范围不同，多尺度反演过程中，通过在大尺度上解的邻域附近来扰动，既保证了有效基因不被丢失又使搜索范围是逐尺度缩小，速度基因值的分辨率提高，这样以来，逐尺度提高了解的精度和分辨率.3井间地震层析成像的发展方向由上可以看出，如何提高图像分辨率和成像效果，仍是今后提高井间地震析成像技术实用化程度所研究的重点.引起图像分辨率偏低、成像效果不好的主要原因是：（1）有限观测角；（2）介质的非均匀性和各向异性；（3）观测数据含有误差；（4）成像方法本身的缺陷（如线性化、正则化和弱散射近似等）.因此，井间地震CT 仍需从反演理论及方法、反演算法，野外施工与数据采集等方面进行研究、改进，进而很好地解决上述几方面的问题，才能提高地震CT 的地质效果.目前，井间地震层析成像技术需要深入研究的主要内容有：（1）含噪声不完全投影数据的3-D 地震层析成像方法；（2）复杂介质井间地震弹性波波场快速数值模拟方法；（3）双相介质井间地震层析成像方法；（4）各向异性介质井间地震层析成像方法；（5）3-D 多波井间地震层析成像方法；（6）矢量（三分量或多分量）井间地震层析成像方法；（7）四维（时间推移）井间地震层析成像方法.参考文献［1］Cerveny V.et al.著，刘福田，等译.地震学中的射线方法［M ］.地质出版社，1986.［2］Van Trier J ，Symes M M.Upwind finite-difference calculation traveltime ［J ］.Geophysics ，1991，56（3）：812～821.［3］周洪波，张关泉.复杂构造区域的初至波走时计算［J ］.地球物理学报，1994，37（4）：515～520.［4］Asakawa E ，Kanawa T.Seismic ray tracing using linear traveltime InterPolation ［J ］.Geophysical Prospecting ，1993，41（l ）：99～111.［5］赵改善，等.基于旅行时线性插值的地震射线追踪算法［J ］.石油物探，1998，37（2）：14～24.・59・第16卷第3期裴正林：井间地震层析成像的现状与进展［6］张霖斌，等.有限差分法射线追踪［J ］.石油地球物理勘探，1993，28（6）：673～677.［7］Qin et al ，Solution of the eikonal equation by a finite difference method ，60th SEG Annual Meeting ，1990.［9］Moser ，T.J.，Shortest path calculation of seismic rays ［J ］.Geophysics ，1991，56（1）：59～67.［8］Vidale J.Finite-difference calculation of traveltlme in three dimension ［J ］.Geophysics ，1990，55（5）：521～526.［9］Moser T J.Shortest path calculation of seismic rays ［J ］.Geophysics ，1991，56（1）：59～67.［10］周兵，赵明阶.最小走时射线追踪层析方法［J ］.物探化探计算技术，1992，14（2）：124～130.［11］Robert L.Coultrip ，High accuracy wavefront tracing traveltlme calculation ［J ］.Geophysics ，1993，58（2）：284～292.［12］刘洪，李幼铭.计算最小走时和射线路径的界面网全局方法［J ］.地球物理学报，1995，38（6）：823～831.［13］范祯祥，等.地震波数值模拟与偏移成像［M ］.河南科学技术出版社，1994.［14］Guust Nolet.冯锐，等译.地震层析技术［M ］.地质出版社，1991.［15］周兵，等.一种新的地震射线层析成像计算方法［J ］.石油物探，1994，33（1）：45～54.［16］吴律，著.层析基础及其在井间地震中的应用［M ］.石油工业出版社，1997.［17］Carrion P.Dual tomography for imaging complex structure ［J ］.Geophysics ，1991，56（9）：1395～1404.［18］张文生，何樵登.约束走时层析成像［J ］.石油地球物理勘探，1997，32（l ）：68～74.［19］殷军，冯锐.井间层析成像的最大滴方法［J ］.地球物理学报，1992，35（2）：234～241.［20］王成礼，孟凡林.井间地震纵横波走时层析成像处理方法［J ］.石油地球物理勘探，1996，31（4）：509～514.［21］Calman et al.Reflection +transmission for crosswell tomographic imaging ［C ］.59th SEG Annual Meeting ，1989.［22］Wu lu et al.Seismic crosshole curved ray reflection plus transmission tomography ［C ］.64th SEG Annul Meeting ，1994.［23］张树林，朱介寿.井间地震和逆VSP 联合层折成像［J ］.石油地球物理勘探，1993，28（5）：577～583.［24］Wu ，R.S.and Toksoz ，M.N.Diffraction tomography and multisource holography applied to seismic imaging ［J ］.Geophysics ，1987，52（l ）：11～25.［25］Lo T W ，Toksoz M N ，Xu S H ，Wu R S.Ultrasonic laboratory tests of geophysical tomographic reconstruction ［J ］.Geophysics ，1988，53（7）：947～956.［26］Pratt R G ，Worthington M H.The application of diffraction tomography to cross-hole seismic data ［J ］.Geophysics ，1988，53（10）：1284～1294.［27］吴律，朱光明.Born 近似层析与Rytov 近似层析对比研究［J ］.石油地球物理勘探，1993，28（1）：1～16.［28］彭成斌，陈隅.衍射CT 技术和多源全息成像技术的比较研究［J ］.地球物理学报，1990，33（2）：154～162.［29］刁顺，等.弹性波衍射CT ［J ］.石油物探，1994，33（3）：47～53.［30］孙豪志，范祯祥.弹性被介质中纵横波速度反演［J ］.石油地球物理勘探，1994，29（6）：678～684.［31］陈湛文，等.弹性介质中密度ρ与拉梅常数λ、µ的衍射层析成像方法研究［J ］.地球物理学报，1995，38（2）：234～242.［32］Pai D M.Crosshole seismic using vertical eigenstates ［J ］.Geophysics ，1990，55（7）：815～820.［33］Dickens T A.Diffraction tomography for crosswell imaging of nearly layered media ［J ］.Geophysics ，1994，59（5）：694～706.［34］黄朕捷，吴如山.垂直非均匀背景多频背向散射层析成像［J ］.地球物理学报，1994，37（1）：87～100.［35］井西利，杨长春.一种非均匀背景的散射层析成像方法研究［J ］.石油物探，1997，36（2）：7～14.［36］Tarantora A.Inversion of seismic reflection data in the acoustic approximation ［J ］.Geophysics ，1984，49（8）：1259～1266.［37］黄联捷，杨文采.声波方程送散射反演的近似方法［J ］.地球物理学报，1991，34（3）：626～634.［38］社剑渊，杨文采.跨孔地震层析成像的级联方法［J ］.地球物理学报，1991，34（6）：771～780.［39］Tura M A C ，Greaves R J ，Beydoun W b.Crosswell seismic reflection /diffraction tomography ：A reservoir characterization applica-tion ［J ］.Geophysics ，1994，59（3）：351～361.［40］Pratt R G ，Goulty N bining wave-equation imaging with traveltime tomography to form high-resolution images from crossholedata ［J ］.Geophysics ，1991，56（2）：208～224.［41］Reiter D T ，Rodi W.Nonlinear waveform tomography applied to crosshole seismic data ［J ］.Geophysics ，1996，61（3）：902～913.［42］王守东，刘家琦.二维声波方程速度反演的一种方法［J ］.地球物理学报，1995，38（6）：833～839.［43］王大川，严忠琼，层析成像中的波形反演［J ］.物探化探计算技术，1996，18（4）：289～303.・69・地球物理学进展2001年［44］Luo Y ，Schuster G T.Wave-equation traveltime inversion ［J ］.Geophysics ，1991，56（5）：645～653.［45］樊计昌，许云.最大摘走时波形层析成像［J ］.石油物探，1995，34（2）：46～52.［46］牛彦良，杨文采.跨孔地震CT 中的逐次线性化方法［J ］.地球物理学报，1995，38（3）：378～386.［47］Bunks ，C ，Saleck ，PM ，Zaleski ，S ，Chavent ，G.Multlscale seismic waveform Inversion ［J ］.Geophysics ，1995，60（5），1457～1473.［48］陈小宏，牟永光.地震波动方程反演的多重网格方法［C ］.中国地球物理学会年刊，地震出版社，1995.［49］钱建良，刘家琦.波动方程反问题的多尺度反演方法［C ］.中国地球物理学会年刊，地震出版社，1995.［50］Zhou C ，Ci W ，Luo Y ，schusteL G ，Hassanzadeh S.Acoustic wave equation trveltime and waveform inversion of crosshole seismicdada ［J ］.Geophysics ，1995，60（3）：765～774.［51］Stoffa P L ，Sen M K.Nonlinear multlparameter optimization using genetic algorithms ：inversion of plane -wave seismograms ［J ］.Geophysics ，1991，56（11）：1794～1810.［52］Sen M K ，Stoffa R L.Nonlinear one -dimentional seismic waveform inversion using simulated annealing ［J ］.Geophysics ，1991，56（10）：1624～1638.［53］Sen M K ，Stoffa P L.Rapid sampling of model space using genetic algorithms ：examples from seismic waveform inversion ［J ］.Geo-phys.J.Int.，1992，108：281～292.［54］Mosegaard K ，Vestergaard R D.A simulated annealing approach to seismic model optimization with spare prior information ［J ］.Geophysical Prospecting ，1991，39（5）：599～611.［55］Pullammanappall ，S K ，Louie J N.Inversion of seismic reflection travetimes using nonlinear optimization scheme ［J ］.Geophysics ，1993，58（10）：1607～1620.［56］石琳柯.地球物理遗传反演方法研究［博士学位论文］［D ］.中国地质大学（北京），1996.［57］Vilis D R ，Ulrych T J.Traveltime tomography using fast simulated annealing ［C ］.65th SEG Annul Meeting ，1995，1055～1057.［58］席道瑛，等.BP 网络的改进及模拟退火神经网络在地学中应用［J ］.物探化探计算技术，1996，18（3）：219～224.［59］杨丽华，孟昭波.一维波动方程小波逐版本反演［J ］.地球物理学报，1995，38（6）：815～822.［60］孟鸿鹰，刘贵忠.小波变换多尺度地震波形反演［J ］.地球物理学报，1999，42（2）：241—248.［61］Zhaobo Meng ，Scales J A.2-D tomography in multi-resolution analysis model space ［C ］.SGE 66th annual meeting ，Expanded Ab-stracts ，1996，1126～1129.［62］裴正林.小波多尺度井间地震层析成像方法研究及其应用［D ］.博士学位论文，北京：中国地质大学（北京），2000.［63］陈国良，等.遗传算法及其应用［M ］.北京：人民邮电出版社，1996.STATUS AND PROGRESS ON CROSSHOLESEISMIC TOMOGRAPHYPEI Zheng-lin（University of Petroleum ，Beijing ，100083，China ）Abstract ：This paper gives a review on crosshole seismic tomography during last decades.And it gives the up-to-date progress on applying the wavelet transform to tomography problem.The paper also gives the prospect of the method.Keywords ：Crosshole seismic tomography ；Wavelet multiscale ；Reserch progress ・79・第16卷第3期裴正林：井间地震层析成像的现状与进展。