（地球探测与信息技术专业论文）从阵列声波数据中反演频散波的相慢度.pdf

石油大学( 华尔) 硕 ：研究生论文 摘要 井下阵列声波测井仪器可以激发并记录在充液井孔中传播的声波， 单极子源激发的斯通利波、伪瑞利波和偶极子源激发的弯曲波都是频散 波。在并孔中，各种频散波的频敌特性在一定程度上依赖周围地层的性 质，它们的相慢度频散曲线是周围地层参数的函数，因此分析声波全波 列中各种频散波的相慢度的频散特性有助于分析周围岩石和地层的性 质，如岩性、密度、各向异性、渗透率、孔隙度、饱和度等，在测井解 释、合成声波波形、设计声波测井仪器等方面都非常有意义。 本文采用p r o n y 方法和阵列声波数据预测处理技术相结合，向前向 后预测得到预测波形，使用快速模拟退火算法求解记录波形和预测波形 残差的全局最小解，反演阵列声波数据中频散波的相慢度频散曲线。本 文将根据频散方程正演得到的频散波的相慢度频散曲线与从阵列声波数 据中反演的结果进行了对比，根据两者的一致性验证了反演结果的真实 性和准确性，并且使用该反演方法处理了连续深度的现场阵列声波波形 数据，得到了满意的结果。另外，本文研究了充液井孔中多极子源激发 的声场，合成了多极子阵列声波全波列。 关键词：频散相慢度p r o n y 方法预测处理快速模拟退火法 石油大学( 华东) 硕士研究生论文 i n v e r s i o no fp h a s es l o w n e s so f d i s p e r s i v ew a y e s f r o m a r r a y a c o u s t i cd 矗t a a b s t r a c t a r r a ya c o u s t i cw e l ll o g g i n gt o o l sc a l le x c i t ea n dr e c o r da c o u s t i cw a v e s p r o p a g a t i n g i naf l u i d f i l l e db o r e h o l e t h e s t o n e l e y w a v ea n d p s e u d o r a y l e i g hw a v ee x c i t e db yam o n o p o l es o u r c e ，t h e f l e x u r a lw a v e e x c i t e db yad i p o l es o u r c e ，a r cd i s p e r s i v ew a v e s i nab o r e h o l e ，t h ep h a s e s l o w n e s sd i s p e r s i o nc u r v e sa r ea f f e c t e d ，t os o m ee x t e n t ，b yt h ep r o p e r t i e so f t h ef o r m a t i o ns u r r o u n d i n gab o r e h o l e s ot h ea n a l y s i so ft h e i rc h a r a c t e ri s h e l p f i l lt od i s c o v e rt h ep r o p e r t i e so f t h es u r r o u n d i n gf o r m a t i o na n dr o c k ，s u c h a sl i t h o l o g y , d e n s i t y , a n i s o t r o p y , p o r o s i t y , p e r m e a b i l i t y , s a t u r a t i o n ，e t c i nt h i s p a p e r , w ei n t e g r a t e d t h e p r o n ym e t h o d a n dt h e p r e d i c t i v e p r o c e s s i n gt e c b m o l o g yo fa r r a ya c o u s t i cw a v e f o r md a t a ，g o tt h ep r e d i c t e d w a v e f o r m s u s i n gt h ef o r w a r d - b a c k w a r dp r e d i c t i o n ，m i n i m i z e dt h ed i f f e r e n c e b e t w e e nt h ep r e d i c t e dw a v e f o r m sa n dt h er e c o r d e do n e su s i n gf a s ts i m u l a t e d a n n e a l i n ga l g o r i t h m ，e s t i m a t e d t h e d i s p e r s i o np h a s e s l o w n e s sc u r v e so f d i s p e r s i v e w a v e sf r o m a r r a y a c o u s t i cw a v e f o r md a t a i nt h i s p a p e r , w e c o m p a r e d t h e s e p h a s ed i s p e r s i o n c u r v e sw i t ho n e sw h i c ha r ec a l c u l a t e d a c c o r d i n gt ot h ed i s p e r s i o ne q u a t i o n t h e i rc o n s i s t e n c yp r o v e dt h et r u e n e s s a n da c c u r a c yo fr e s u l t s ，w eu s e dt h i sm e t h o dt o p r o c e s st h e f i e l d a r r a y a c o u s t i cw a v e f o r md a t ai nc o n t i n u a ld e p t hp o i n t sa n dg o ts a t i s f a c t o r yr e s u l t s i na d d i t i o n ，w es t u d i e dt h ea c o u s t i cf i e l de x c i t e db yt h em u l t i p o l es o u r c ei na f l u i d f i l l e db o r e h o l ea n ds y n t h e s i z e da r r a r ya c o u s t i cw a v e f o r m s k e yw o r d ：d i s p e r s i o np h a s e s l o w n e s s p r o n y m e t h o d p r e d i c t i v e p r o c e s s i n g f a s ts i m u l a t e d a n n e a l i n gm e t h o d 石油大学( 华东) 硕士研究生论文 独创性声明 本人声明所呈交的论文是我个人在导师指导下进行的研究工作及取 得的研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论 文中不包含其他人已经发表或撰写过的研究成果，也_ i 包含为获得石油 大学或其它教育机构的学位或证书而使用过的材料。与我同工作的同 志对本研究所做的任何贡献均已在论文中作了明确的说明并表示了谢 意。 签名：兰篷塞：国 又伊奶年岁月日 关于论文使用授权的说明 本人完全了解石油大学有关保留、使用学位论文的规定，即：学校 有权保留送交论文的复印件及电子版，允许论文被查阅和借阅；学校可 以公布论文的全部或部分内容，可以采用影印、缩印或其他复制手段保 存论文。 ( 保密论文在解密后应遵守此规定) 学生签名： 导师签名： 2 0 0 3 仨5 2 0 0 3 芷5 月6 日一 月6 口一 4 i 油人学( 华东) 硕士研究生论文第1 章引言 第1 章引言 1 1 论文研究目的及意义 1 1 1 目的 声波测井仪器可以激发并记录在充液井孔中传播的声波，记录的波 形中包括不频散的纵波和横波，以及单极声波测井仪激发的频散的斯通 利波、伪瑞利波，偶极声波测井仪激发的频散的弯曲波。频散是指波在 不同的频率点处有不同的相慢度和群慢度( 慢度即速度的倒数) 。各种频 散波的频散特性在一定程度上依赖周围岩石的性质，具体地说地层的性 质不仅影响声波的幅度而且影响它的慢度，并且这种影响在不同频率是 不同的。因此需要使用信号处理技术从阵列声波波形数据中估计各种频 散波在不同频率的慢度，从阵列声波波形数据中提取频散波的慢度。 1 1 2 意义 在井孔中，各种频散波的频散特性在一定程度上依赖周围地层的性 质，它们的相慢度频散曲线是周围地层参数的函数，如各种波在地层中 的传播速度、地层的密度、井孔中泥浆的密度。因此反演频散波的相慢 度频散曲线对于从波形资料中准确获得这些参数有很重要的作用，特别 在从偶极子波形资料中获得地层横波慢度时，为了准确提取横波慢度， 必须考虑弯曲波的频散，进行频散校正，这时必须用到弯曲波的相慢度 频散曲线。另外只有准确的得到了波速、密度等地层参数，才有可能根 据各种模型反演地层的渗透率、孔隙度、各种饱和度，如利用威利和雷 曼公式根据纵波的慢度求孔隙度，根据斯通利波的慢度求渗透率。因此 分析声波全波列中各种波的频散特性有助于分析周围岩石和地层的性 质，如岩性、密度、各向异性等，在测井解释、合成声波波形数据、设 石油大学( 华东) 硕士研究生论文第1 章引言 计声波测井仪器、储层评价等方面都具有重要意义。 1 2 国内外研究现状 1 2 1 概论 阵列声波测井的波形数据通常是由多种波动模式构成，如单极全波 列测井的波形数据是由纵波、横波、伪瑞利波、斯通利波构成，偶极横 波测井的波形数据主要由弯曲波构成。如果这些波在时域是分离的，就 可以使用传统的波形处理方法处理，如各种滤波、n 次根、相关的方法 得到纵波、横波、斯通利波的慢度。然而在很多情况下，这些波动模式 在时域和频域相互重叠，如在软地层的单极全波列测井中，横波和斯通 利波相互重叠，在各向异性地层的偶极横波测井中，慢横波和快横波相 互重叠。这时传统的处理方法就无能为力了，而且传统的处理方法处理 时没有考虑频散波如弯曲波和斯通利波的频散特性，因此如何从阵列声 波波形数据中将这些波分离并且反演频散波的慢度成为一项很有意义的 _ _ r + - 一。 1 2 2 现状 目前研究从阵列声波波形数据中反演频散波的慢度的国内外文献比 较少，以下介绍一下从事这方面研究的理论和这些理论的主要思想。 乔文孝和阎树汶( 1 9 9 5 ) 使用相位展开法提取斯通利波的慢度【1 】， 这种方法只需两个通道的信号即可。相位展开法首先选取两个通道信号 中包含频散波的波形部分开窗，然后用自适应数值积分的方法对窗中的 信号分别进行相位展开，最后计算对应每一频率的相位差、延迟时间及 相速度。这种方法实用，但没有考虑波的衰减，没有充分利用阵列声波 2 石油大学( 华东) 硕士研究生论文第1 章引言 波形数据，达到降低噪声干扰，提高准确性的目的。h s u 和b a g g e r o e r ( 1 9 8 6 ) 使用最大似然的方法从阵列声波波形数据中估计各种频散波在 不同频率的慢度 2 。这种方法得到的结果准确，但计算非常耗时。h s u 和e s m e r s o “1 9 9 2 ) 使用参数估计的方法估计各种频散波的相慢度和群慢 度 3 1 。他们近似的认为各种波的波数是频率的线性函数，预测某一接收 探头接收到的波形，然后最小化预测得到的和真实记录的波形，从而估 计各种波的相慢度和群慢度。以上两种方法都可以较准确的估计各种频 散波的相慢度，但它们都没有考虑各种频散波的衰减，都涉及到复杂的 矩阵的各种计算，求解算法较难实现，计算量大。 l a n g 等( 1 9 8 7 ) 使用著名的指数模型的p r o n y 理论估训各种频散波 的慢度【4 】，它又称为向前一向后线性预测方法该方法根据井孔中各种波 的传播规律得到p r o n y 指数模型，然后使用p r o n y 方法向前一向后预测 某一接收探头接收到的波形，最后将预测得到的波形和实际记录的波形 进行比较，从而估计各种频散波的在不同频率的慢度，l a n g 等并没有考 虑各种频散波的衰减。e l l e f s e n ( 1 9 8 9 ) 也使用扩展的p r o n y 方法不仅得 到了阵列声波波形数据中各种频散波在不同频率的相慢度而且得到了它 们的衰减【5 】。 图1 1 是l a n g ( 1 9 8 7 ) 使用p r o n y 方法处理合成的单极阵列声波全 波列波形数据。 石油大学( 华东) 硕士研究生论文第1 章引言 00 5 11 522 53 t i m e ( m s ) 图1 1 合成的单极阵列声波全波列波形 图1 1 是在变源距，等间距的情况下合成的4 0 道单极全波列波形， 从波形中可以看到幅度微弱的纵波，然后是幅度较强的横波，最后是幅 度最强的斯通利波。 4 打油大学( 华东) 硕士研究生论文第1 章引言 要 差 p o w e rs e a i e -0d 日 一2 0d b 一4 0d b 一6 0d b - 一8 0d b f r e q u e n c y ( h z ) s i g n a 匕r r o r 图1 2 p r o n y 方法的处理结果 图1 2 a 是各种波的相慢度频散瞳线，从图中可以看出慢度随频率的 变化规律。慢度最小的是纵波，它几乎不频散。然后是三种模式波，它 们都有很明显的频散，存在截至频率。慢度最大的是斯通利波，它也有 频教，但没有截至频率。图1 2 b 是各种波的波数和频率的变化关系。图 1 2 c 是信号和残差随频率变化的曲线。 一e 、3 mmc言一 一#、st，)们c言is lde丁ce，|口mni|meloz 俐 埘 石油人学( 华东) 硕士研究生论文第1 章引言 f r e q u e n c y ( k h z ) 图1 3 理论计算的相慢度频散曲线 图l - 3 是理论计算的相侵度频散曲线。对比图13 和圈1 2 a 可以看 出理论计算的相慢度频散曲线和p r o n y 方法得到的相慢度频散曲线是一 致的，这说明p r o n y 方法得到的结果是可信的。 1 2 3 小结 由于采用p r o n y 方法从阵列声波波形数据中反演频散波的相慢度频 散曲线这项波形处理技术可以不受不同波动模式在时域重叠的影响，而 且可以反演出频散波在各个频率点的慢度，从而分析各种频散波的频散 特性。p r o n y 方法得到的结果准确可靠，不但可以反演各种波的慢度也可 以反演衰减，而且它的计算效率很高。目前这项波形处理技术在国内还 没有具体实现的文献，国外虽然有研究，但是没有将这项技术应用于处 理现场的阵列声波波形资料，因此对这项技术的研究和探索是非常有意 义和价值的。 1 3 论文的研究内容 本论文实现了将阵列声波波形数据预测处理技术和p r o n y 方法结合， 使用快速模拟退火算法从阵列声波波形数据中反演各种频散波的相慢 6 一、slo詈ls 富、搴i，蕾譬i-t耍 石油大学( 华东) 硕士研究生论文第1 章引言 度，得到了合成阵列声波波形中的弯曲波、斯通利波、伪瑞利波和现场 阵列声波波形中的弯曲波、斯通利波、套管波的相慢度频散曲线，并将 反演结果与根据频散方程m 演得到的相慢度频散曲线进行了对比，两者 的一致性验证了该方法的准确性，得到了满意的结果。另外本沧文研究 了充液井孔中多极子源激发的声场，利用实轴积分的方法合成了多极予 源激发的阵列声波波形。 7 石油大学( 华东) 硕士研究生论文 第2 章p r o n y 方法和预测处理技术 第2 章p r o n y 方法和预测处理技术 2 1 p r o n y 方法及其在阵列声波测井中的应用 2 1 1 p r o n y 方法的理论基础 p r o n y 方法是用一组指数项的线性组合来拟合等间距采样数据的方 法，可以从中分析出信号的幅值、相位、阻尼因子、频率等信息。每个 指数项包含幅度和指数因子两个部分，对于个n 阶的模型，就含有知 个待求因子，通常测量点数取为n 2 n ，并用最小二乘的估算方法，减少 噪声的影响。 p r o n y 方法用, 个任意幅值、相位、阻尼因子、频率的指数模型来拟 合被测数据，即假定模型具有如下形式： y 。( 七) = b i z i ( 七= o 1 ，n 一1 ) ( 2 1 ) 其中，旦，互为复数，包含了确定信号频率和阻尼因子信息的量，其一般 形式可写为： f e = a ie x p ( j o j ) 【z i = e x p ( a + j 2 z f ) a t 式中4 一振幅 醴一相位 口i 一衰减系数 z 一振荡频率 ( 2 2 ) 加一时间采样间隔 为了求出4 、q 、吒、z ，使拟合数据和被测数据的误差平方和 8 石油大学( 华尔) 硕p 研究生论文 第2 章p r o n y 方法年预测处理技术 占= i y ( 女) 一y ( 州 ( 2 3 ) 最小即可求得，但这将成为一个复杂的非线性最小_ _ ：二乘问题。p r o n y 方法 的关键是认识到y 。( ) 是一常系数线性差分方程的解，即： y ( 尼) = 一n 。y ( 一m ) ( 2 4 ) m = l 其中t ；。是如下方程的根。 兀( z - z d = q z n - i , = l = li = 0 ( 2 5 ) 因此只需求得系数a 。，即可通过解多项式的根得到指数参量。设 则： 令 于是 e ( k ) = y ( k ) - y ( 七) ( 0 k n 一1 ) _ y ( ) = 一d 。y ( k 一所) + e ( k - m ) s ( ) = e ( 后一m ) y ( 七) = 一y ( k - m ) + c ( k ) ( 2 6 ) ( 2 7 ) ( 2 8 ) ( 2 9 ) n - l 要求系数，只要使is ( 七) 1 2 最小，这样就可以用最小二：乘估计得 k = n 到。知道后，就可以求得多项式的根。由式( 2 1 ) 可以得到以b i 为未 9 互塑盔堂哆东) 硕士研究生论文第2 章p r o n y 方法和预测处理技术 知数的线性方程组 z ? z 乏 z ! 霉 z 掣。1z y 掣1 且 垦 ： b ， y ( o ) y ( 1 ) y ( n n ( 2 1 0 ) 使用最小二乘可以得到e 。得到旦，z 。后，就可以确定信号的振幅、相位、 阻尼因子和频率。 随着当今计算机技术的发展，p r o n y 方法得到了越来越广泛的重视和 应用。在数字信号处理中，对于由指数线性组合而成的信号，p r o n y 方法 相对于常用的信号处理技术如傅旱叶变换有如下优点： ( 1 ) p r o n y 方法可以直接求得信号的幅值、相位、阻尼比和频率，而无 需从频率响应来求，所以计算量大为减少。 ( 2 ) p o m y 方法采用最小二乘法意义上的拟台，有利于消除噪声的影响。 2 1 2p r o n y 方法用于阵列声波测井信号 将p r o n y 方法应用于阵列声波波形数据的处理，首先应当根据声波 在井孔中的传播规律建立阵列声波测井模型来拟合记录的波形数据，然 后根据p r o n y 方法估计阵列声波测井信号的慢度等参数。 2 1 2 1 阵列声波测井模型的建立 p r o n y 方法使用指数模型估计在各个频率上的空间波数，它又被称为 向前和向后线性预测。指数模型是从井孔中声波传播的物理规律得到的， 它是声波信号处理的基础。 首先假设与接收阵列在同一深度的地层性质在井轴方向上不变，径 向上可以任意改变。根据t s a n g 和r a d e r ( 1 9 7 9 ) 【6 ，通过接收阵列的声 波在井孔的轴向上的传播可以式( 2 1 1 ) 描述。 1 0 堑迪查堂! 坐东) 硕士研究生论文第2 章p r o w 方法和预测处理技术 w ( z ，f ) = j x ( ) 爿( 七，( d ) e - i c o t e j k z d k d c o ( 2 1 1 ) 其中似z ，t ) 表示距声源距离z ，时刻t 的井轴上的声压变化。a ( k ，) 表示 地层对应的频率一波数响应，k 代表垂直方向的波数，0 9 表示角频率。 x ( c o ) 表示发射脉冲的傅里叶变换。 根据式( 2 1 1 ) 可以得到w ( z ，t ) 在频域中的表达式 w ( z ，c o ) = x 沏) f a ( k ，c o ) e m d ( 2 1 2 ) 上式对波数k 的积分包括声波全波列中各种波的贡献。根据t s a n g 和 r a d e r ( 1 9 7 9 ) 州，k u r k j i a n ( 1 9 8 5 ) f 7 】对波数j i 的积分是a ( k ，c o ) 的孤 立点的贡献之和。这些孤立点有两种：极点和支点。极点的贡献与模式 波相联系，支点与体波相联系，因此式( 2 1 2 ) 也可以表示为： 肜( z ，c o ) = 2 丌j ( 极点的留数) + 且支割线) ( 2 1 3 ) 下面分别讨论极点和支点在式( 2 1 3 ) 中的贡献。 对于位于k = ( 国) 的极点，极点的留数是( 1 2 丌j ) ( c o ) e x p j k p ( ) z 的形式。将它代入( 2 1 3 ) 式得到 矿( z ，国) = a p ( c o ) e x p j k ，( ) z 】+ 且支割线) ( 2 1 4 ) p 。l 在声波测井应用中，极点在式( 2 1 4 ) 中的贡献可以通过以下方法 简化。由于接收探头一般放在声源的近声场之外，因此实的k 主导着虚 的七。这可以通过( 2 1 4 ) 式看出，因为实的。不会随z 的增加而衰减， 而虚的k 。随z 的增加呈指数衰减。因此对于给定的c o ，存在有限个实的 石油大学( 华东) 硕士研究生论文第2 章p r o n y 方法和预测处理技术 极点，无限个虚的极点，故在声波测井中，对测量到的信号的贡献主要 是有限个实极点，剩余的可以忽略。式( 2 1 4 ) 又可以写为 矿( z ，c o ) = 警a p ( ) e x p j k ，) z + r 支割线) ( 2 1 5 ) 其中p ( c o ) 是实的极点七。的个数。 现在考虑支点对接收到的信号的贡献。在数学上割线对连续谱的积 分等于支点和附近的极点的贡献之和。在源声场中，如果附近没有极点， 支点对积分的贡献起主导作用( k u r k j i a n ，1 9 8 5 ) 【7 1 。在声波测井应用中， 在发射探头激发的某些频带中，事实上在支点的附近有极点，这就是说 在这些频带内支割线的贡献可以认为是极点的贡献，而在这些频带外， 支割线的贡献相对来说很小可以忽略。 将支割线的贡献认为是极点的贡献，就得到了空间一频率域的模型， 它是声波信号处理的基础。 w ( z ，c o ) * 口p ( c o ) e x p j k p ( 彩) z ( 2 1 6 ) p = l 现在的问题就是估计t 。和a 。，去拟合记录的阵列波形。 2 1 。2 2 阵列声波测井信号的参数估计 p r o n y 方法( 即向前一向后预测) 的实质是当用指数信号的和去拟合 波形数据很困难时( 实际上它是一个非线性最小二乘问题) ，p r o n y 方法 将它转化为用记录的波形数据的线性组合去预测空问某点的波形数据， 拟合该点记录的波形数据，它需解决的问题就是确定用于预测的各记录 波形的系数，实际上这是一个线性最小二乘问题。 根据l a n g ( 1 9 8 7 ) 【4 等，对于给定频率，拟合的波形和记录的波形 的残差是： 1 2 钉油大学( 华东) 硕士婿究生论文第2 章p r o n y 方法艰l 预测处理技术 m l p埘pip 以哟+ “m n ) 6 ( ”) f 2 + f ( m ) + + ) 6 ) f 2 ( 2 1 7 ) m = p n = lm = on = l 其中“卅) = w ( z 。，) ，m = 0 1 ，m 一1 ，w ( m ) 是w ( m ) 的共轭，6 ( h ) 是待定 系数，式( 2 1 7 ) 中第一项表示向前预测的残差，第二项表示向后预测 的残差。接下来的问题就是如何使残差最小，这足一。个线性最小二乘问 题。在确定了6 ( n ) 后，最小化式( 2 1 8 ) 就可以得到k p 的估计值，然后 根据屯= k p 得到慢度。 b ( k ) 1 2 = 1 1 + b ( n ) e 肼1 2 ( 2 1 8 ) 向前一向后预测方法有较好的稳定性，而且不要求各种波在时域里 是彼此分离的，但这里并没有考虑波数的虚部产生的衰减，即波数是实 的。另外这种方法在估计过程中会产生虚假的波数值，所以有必要对得 到的波数估计值进行编辑，舍去虚假值。可以采用能量编辑的方法，舍 去那些对应能量l a 。1 2 很小的波数。 2 2 阵列声波数据的预测处理技术 2 2 。1 简介 阵列声波测井数据通常是由多种波动模式构成，如单极全波列波形 数据是由纵波、横波、伪瑞利波、斯通利波构成。由于这些波的慢度、 频率成分是可以区分的，因此可以使用传统的波形处理方法处理，如各 种滤波、n 次根、相关的方法。然而在很多情况下，许多波动模式在时 域和频域相互重叠，如在套管井的偶极横波测并中，偶极子波和套管波 相互重叠，在各向异性地层的偶极横波测井中，慢横波和快横波相互重 叠，使用普通的方法很难将它们区分开。而预测的处理技术可以用来处 1 3 石油大学( 华东) 硕士研究生论文 第2 章p r o n y 方法和预测处理技术 理包含多种波动模式的声波阵列数据。它使用时域里的预测理论，通过 组合其它接收探头的波形来拟合另外一个接收探头的波形，然后最小化 记录的波形与预测的波形之差，从而达到估计这些波的慢度的目的。最 重要的是这种最优化过程是直接在时域里实现的，因此很有效地减少了 噪音的干扰，增强了准确性。 2 2 2 理论和方法 阵列声波数据的预测处理技术( t a n g ，1 9 9 7 ) 【8 】的理论基础是线性 预测理论。假设个接收探头接收到的波形由p 种波动模式构成，根据 线性预测理论，第一个接收探头的接收到的波形的频谱可以用其它几个 接收探头的频谱的线性组合来表示，即： p r v ( c o ) = 一d 。( 国) 。( 国) ，o = 1 ，) ( 2 1 9 ) = 1 系数口与特征公式( 2 2 0 ) 相联系。 z + 口l z 一1 + t + 口p 一】z + a p = o ，( z = e 删) ( 2 2 0 ) 其中d 是接收探头的间距，s 表示慢度。式( 2 1 9 ) 和式( 2 2 0 ) 是p r o n y 方法的基础。它的总体思想是根据式( 2 1 9 ) 得到系数a ，代入式( 2 2 0 ) ， 通过确定方程的根来估计慢度。预测处理技术将系数a 用式( 2 2 0 ) 的 根来表示。 口1 = 一z n ：= = 。z ， ( k ( ，) n ，= 一z ；： k ，m 1 a p = ( 一1 ) 9z 1 = 2 zp 1 4 ( 2 2 1 ) 石油大学( 华东) 硕士研究生论文 第2 章p r o n y 方法和预测处理技术 将式( 2 | 2 1 ) 代入式( 2 1 9 ) ，这样被预测的波形就成为p 种波的慢度的 函数，如果忽略慢度的频散，将式( 2 1 9 ) 傅罩叶反变换到时域，则得 到预测的接收探头n 接收到的波形。 睨( f ) = + ( f 文d ) 一。：( f ( + j ) d ) ( ) ( 2 2 2 ) + + ，( f ( & + s + ) d ) + 一( 一1 ) ，( f d ) k = l ，m = l k = l ( 0 。国”代表角频率的虚部，c o 代表角频率的实部。具体的积分 方法可以用辛普生二维积分公式。引入复角频率和非负的“的目的是为 了确保沿着积分路径没有地层响应函数爿( 屯，) 的奇点。如果0 9 ”为零， 常见的波动模式的极点就会落在实t 轴上。 在进行二维积分时，先对t 积分，然后对进行积分，并且采用下 列规则确定积分的采样间隔触：和国，积分的范围( t ) 一和( 。) 。 ( 1 ) 确定t 的积分采样间隔屯：t = m i n ( 国”3 也，0 2 z ，0 2 a ) 。 ( 2 ) 确定t ：的积分范围( 也) 一：( 也) 。= m a x ( 5 1 a ，1 5 ：o 7 v ，) ，那么在t 积 分时的采样次数m 是：帆= 訾+ 1 。 ( 3 ) 确定+ 的积分采样间隔a o o ：a c o 。= m i n ( o 2 0 f , 0 2 m ) 一，0 2 v ，d ) 。 其中( f ) 一是合成波形的持续时间。 ( 4 ) 确定的积分范围( + ) 。：( c o ) 一= 5 。其中0 3 。是中心角频率。 那么在珊积分时的采样次数心：。= 号手+ l 。 ( 5 ) 确定角频率的虚部”：0 9 ”= 2 ( t ) 一。 石油大学( 华东) 硕士研究生论文第3 章充液井孔中的多极子声波 在数值计算中，使用b l a c k m a n 窗函数的_ 二阶导数作为源函数 ( k u r k j i a n ，1 9 8 3 ) 1 1 1 , b l a c k m a n 窗函数和它的二阶导数为： 们) ：接和s ( 了2 r r n t ) o c 。 c r 或 ，。( f i t ) 和k 。( 唧) ，j 。( p r ) 和 口2 0k ( ) c 。 c r 或a2 判断是否接受= 接受或舍弃”的迭代，并逐步衰减t 值，算法终止时的当 前解即为所得近似最优解。退火过程由冷却进度表控制，包括控制参数 的初值t 及其衰减因子f 、每个t 值时的迭代次数( 称为一个m a p k o b 链的长度) l 和停止条件鼠 模拟退火算法的主要优点是：不用求目标函数的偏导数及解大型矩 阵方程组，就能找到一个全局最优解，而且易于加入约束条件，编写程 序简单。 4 。1 。2 基本原理 模拟退火算法( s a ，s i m u l a t e d a n n e a l i n g ) 是一种启发式的蒙特卡罗 ( m o n t ec a r l o ) 方法。这种算法是在给定的模型空间内搜索目标函数达 到全局极小值的最优模型，它已用于各种最优化问题的计算。常规的模 3 6 石油火学( 华东) 硕士研究生论文第4 章反演频散波慢度的程序实现 拟退火算法即k i r k p a t r i c k 等( 1 9 8 3 ) 提出的m e t r o p o l i s 算法，包括以 卜- 步骤： ( 1 ) 给定模型每一参数的变化范围，在这个范围内随机选择一个初始模 型，并计算相应的目标函数值e ( m 。) 。 ( 2 ) 对当前模型进行扰动产生一+ 个新模型m ，计算相应的目标函数 值e ( m ) ，得到e = e ( m ) 一e ( m 。) 。 ( 3 ) 若业 0 ，则新模型m 按概率 p = e x p ( 一a e i t ) 进行接受，丁为温度。当模型被接受时，置 m o = m ，e ( m o ) = e ( m ) 。 ( 4 ) 在温度，下，重复一定次数的扰动和接受过程，即重复步骤( 2 ) 和( 3 ) 。 ( 5 ) 缓慢地降低温度r 。 ( 6 ) 重复步骤( 2 ) 、( 5 ) ，直到满足收敛条件为止。 以上的算法实际分两步交替进行计算：随机扰动产生新模型并计算 目标函数的变化；决定新模型是否被接受。由于算法是在高温条件开始 进行的，因此使e 增大的模型可能被接受，这样就能舍去局部极小值。 通过缓慢地降低温度，算法能收敛到全局最优点。 非常快速的模拟退火算法( v f s a ) 具有较高的计算效率其主要原因 在于用依赖于温度的似c a u c h y 分布产生新的扰动模型。有关v f s a 的模 型扰动、降温方式如下： 模型扰动 模拟退火算法中新模型的产生是对当前模型进行扰动得到的。i n 曲e l ( 1 9 8 9 ) 在其提出的非常快速的模拟算法中1 4 】，采用依赖于温度的似 3 7 石油大学( 华东) 硕士研究生论文第4 章反演频散波慢度的程序实现 c a u c h y 分布产生新模型，即： ；= m 。+ y ，( 马一4 ) y i = t s g n ( u o 5 ) 【( 1 + 1 t ) 1 2 “。1 【一1 式中：m 。为当前模型中的第i 个变量；u 为( o ，1 ) 均匀分布的随机数； 4 ，茸 为的取值范围，且要求扰动后的m ： 4 ，b i ；s g n ( x ) 为符号函数。采用 依赖于温度的似c a u c h y 分布产生的新模型的优点是：在高温下进行大范 围的搜索，在低温时搜索仅在当前模型附近进行，而且由于似c a u c h y 分 布有一平坦的“尾巴”，使其易于迅速跳出局部极值。这一改进加快了模 拟退火方法的收敛速度。 降温方式 i n g b e r ( 1 9 8 9 ) 1 1 4 给出的非常快速模拟退火方法的降温方式 r ( ) = t oe x p ( 一c k “m ) 式中：t o 为初始温度；k 为迭代次数；c 为给定 常数；n 为待反演的参数的个数。上式可改写成r ( k ) = 瓦口“。 3 8 石油大学( 华东) 硕士研究生论文第4 章反演频散波慢度的程序实现 4 1 3 程序流程图 图4 1 说明了计算机实现快速模拟退火算法的基本步骤。 快速模拟退火算法程序流程示意 注释：n s u b it e r a i o r s 是子步迭代最大次数 开始 i 设置初值模型m o = x o j 计算目标函数值f ( m o ) 0 对当前模型m 0 进行随机扰动得到m 1bj- 计算目标函数值f ( m 1 ) 上 d f = v ( m i ) 一f ( m o ) 喇上匠一， 设m o = m 1 f ( m o ) = f ( m 1 ) 一一一：f 。反丽平面磊l 磊毒军忑忑i ；i i 而而一 一一工堡坚堕坚璺苎! 尘! ! ! 竺! 竺竺型翌l 一一一 喇0 按照快速模拟退火算法的降温方式降温 一一一一一一五幕磊：西藏一一 一一苎竺竺兰! ! 一一一一一 划0 保存f 全局摄小值和对应的参数值 0 终止 图4 1 快速模拟退火算法程序流程图 3 9 石油大学( 华东) 硕士研究生论文第4 章反演频散波慢度的程序实现 4 2p r o n y 方法和预测处理技术结合及实现 根据第二章中的p r o n y 方法和预测处理理论，l a n g 等( 1 9 8 7 ) 使用 指数模型的p r o n y 理论从阵列声波波形数据中反演各种频散波的慢度【4 】， 这种方法的实质是在频率域中，使用最小二乘方法最小化根据向前向后 预测得到的预测波形和记录波形的残差。而t a n g ( 1 9 9 7 ) 阵列声波波形 数据处理理论也是建立在p r o n y 理论的基础上 8 】，不仅如此，预测理论还 得到了特征方程系数的解析表达式，它是以特征方程的根表达的，而且 得到了时域中向前向后预测而得到的预测波形的表达式。即： 睨( ) = 吸。( t - l - _ s k d ) + 艺呒。，o 艺( 艮+ o + s 。) d ) + 一( 一1 ) 呒+ ，( t + - z s i d ) ( k , j , j m = m 1 ) 2 l 式中d 是接收探头的间距，p 表示波形由p 种波动模式构成，s 表示慢度， 矿。是预测波形。 本文将p r o n y 方法和预测处理技术结合，从阵列声波波形数据中反 演频散波的慢度。首先对时域中的预测波形做傅里叶变换，得到频率域 中预测波形的表达式即： 哟巩阳) 抄叫( 珊) 啪4 ) = + ，( 国) e “叫+ ：( 珊) p 张“ 2 2 12 1 ( 4 2 ) 口 三“。未，( s 1 + s y + s m ) d 兰叫 + 吸。， ) 艺p o 堀 + 一( 一1 ) 呒。，( ) p “r ，然后在每个频率点上使用快速模拟退火算法最小化记录波形与预测波 形的残差，反演各种波的慢度，最后将频散波在各个频率点上得到的慢 4 m 0 + 唧 ，粥 o 蛆吸 ，h 石油大学( 华东) 硕士研究生论文 第4 章反演频散波慢度的程序实现 度连接起来就得到各种频散波的相慢度频散曲线。使用该方法可以避免 在l a n g 等( 1 9 8 7 ) 的方法中为了消除虚假结果而对结果进行的繁琐的能 量编辑1 4j ，另外本文