基于偶极弯曲波频散的横波慢度径向分布反演

马明明，陈 浩，何 晓，王秀明

声场声信息国家重点实验室 中国科学院声学研究所，北京 100190

1 引 言

地层横波慢度信息对石油勘探开发有着重要的意义，它与地层岩石性质密切相关.井壁周围地层性质常常因为井壁应力集中、钻井机械损坏和泥浆入侵等外力因素发生变化，进而导致井壁周围地层出现横波慢度的径向变化带.反之，如果已知径向变化带横波慢度分布情况，就可以了解井壁周围地层性质的改变.对于慢地层（地层横波慢度大于井孔流体纵波慢度），传统的单极子仪器探测不到横波首波，原因是井外地层的横波慢度大于井孔流体纵波慢度时，全波列中无法探测到横波首波信号.因此，人们通常采用偶极弯曲波来反演地层横波的信息［1－2］.

1996—2006年Sinha等发展了基于微扰法［3］和BG理论［4］的反演算法［5－8］，来反演地层横波慢度的径向分布，在处理实际数据时，此方法不能同时满足高分辨率和高精度需求.2000年，Braunisch等又提出一种利用多频、多模式波的信息联合反演地层横波慢度的方法［9］，但此方法在慢地层应用受限.2010年，Tang等建立了参量化的、高频约束条件下的反演方程［10］，方程假定地层横波慢度是径向半径r的指数函数，并且在反演的过程中对弯曲波进行高频约束，但高频弯曲波的激发强度较小，易受噪声等的干扰，为反演结果的准确性带来了挑战.2011—2012年，Yang等建立了套管井模型下的参量化的反演方程［11－12］，他们也假定地层横波慢度是径向半径r的指数函数，指数函数的变化形式由3个参数控制，反演地层横波慢度的径向分布，转换为对这3个参数的反演，但是反演计算方法的收敛性受初始值影响严重.

另外，国内的陈雪莲和王瑞甲研究了TI地层模式波对井孔流体参数和地层参数的灵敏度，分析了各模式波的径向探测深度，其中弯曲波探测深度约为4个井孔半径［13］.这也说明可通过弯曲波来获取井壁周围地层的横波慢度，进而了解地层的性质.

前人的工作都说明了，利用弯曲波频散可以求解井壁周围地层横波慢度的径向分布，进而了解井壁周围的地层性质.但是，前人对所采用的反演方法的特性分析较少，并很少涉及如何提高计算效率.本文的工作则主要侧重于这两个方面.本文在前人的工作基础上，首先计算了地层横波慢度径向均匀（以下简称均匀模型）和非均匀情况下（以下简称非均匀模型，非均匀模型的产生是由于均匀模型的横波慢度在径向上发生了变化）的时域全波列波形.其次利用频域方法［14］提取了波形的弯曲波频散曲线，并对弯曲波频散曲线进行了拟合校正［15］.再次，讨论了如何结合微扰法与BG理论，基于非均匀模型和均匀模型之间弯曲波的差异，进行地层横波慢度径向分布的反演计算.最后，我们分析了如何在反演的迭代过程中提高计算效率，并且用多种计算实例，验证了所采用的反演算法的鲁棒性，为它的现场应用打下基础.

2 偶极井孔声场理论

在三维柱坐标系（r，θ，z）下，以偶极源的中心为原点，并使z轴与井轴方向一致.根据文献［16］，井孔内位移势函数φf在频率波数内的稳态声场表达式如下：

井外地层声场纵波位移势函数，SH横波位移势函数和SV横波位移势函数在频率波数内的稳态声场表达式如式：

其中，n＝1代表偶极源，ra为偶极距，r为径向半径，εn＝2，ω为角频率，In和Kn为n阶虚宗量贝塞尔函数，A′n、An、Bn、Cn、Dn、En和Fn是待定系数，由井壁处的边界条件确定分别为井孔流体纵波、井外地层纵波和井外地层横波的径向波数，sf、sp和ss分别为井孔流体纵波、井外地层纵波和井外地层横波慢度，k是轴向波数.

在相邻2个层的界面上声场应满足边界条件.如果在半径rj（下标j＝0，…，N表示径向上的层的序号，N表示径向上的层数）处的界面两侧都是固体，在界面上胶结良好，则界面上的位移和应力分量应该连续，即

式中，ur为法向位移分量，uθ和uz为切向位移分量，σrr为法向应力分量，σrz和σrθ为切向应力分量.

在井壁r0处，其界面的内侧是液体，外侧是固体，此时边界条件变成法向位移ur和法向应力σrr连续，切向应力σrz和σrθ为零，切向位移分量uθ和uz无需考虑，即

层与层之间的声场由汤姆森—哈斯克（Thomson－Hanskell）传播矩阵［17］来连接，反复利用各层介质的传递矩阵和边界条件，就可以得求的径向分层地层的声场.对于偶极声源，测量的是径向位移，是瞬态声场，需要考虑声源的影响，其在空间－时间域的具体表达式为［16］

2.1 时域波形

计算分析时，本文给出了两个模型，模型1的井外是无限大均匀地层，模型2的井外地层横波慢度径向分层，它是由模型1的地层横波慢度径向扰动所致.两个模型的井孔半径r0＝0.1m，井孔流体密度ρf＝1000kg／m3，井孔流体纵波慢度sf＝666.7μs／m.模型1的地层参数见表1，模型2中，地层横波慢度径向是非均匀的，导致剪切模型模量径向非均匀，进而引起纵波慢度径向非均匀，然而纵波慢度对弯曲波频散影响很小［3］，所以为了提高计算效率，计算过程中假设模型2的地层纵波慢度与模型1相同，此外，地层密度对弯曲波频散影响也很小，在计算时也同样假设模型2的地层密度与模型1相同，模型2中地层横波慢度ss的径向分布数值见表2.图1为模型2的示意图，其中径向半径b和井壁r0之间为横波慢度的变化带，b＝0.3m.

表1 地层参数Table 1 Parameters of formation

表2 变化带横波慢度值Table 2 The values of shear wave slowness of altered zone

图1 井孔与地层模型Fig.1 The model of borehole and formation

文中计算用到的声源函数为F（ω）＝声源中心频率f0＝3kHz，声源带宽fb＝0.4f0，声源与第一个接收器之间的距离d＝3.0m，共八个接收器，相邻接收器之间的距离dz＝0.1524m.

图2中黑色和红色曲线分别表示模型1和模型2的偶极全波波列，从图2可以看出模型2的全波波列明显延后于模型1的全波波列，这是由于模型2的井壁周围地层横波慢度增加所致.

2.2 频散分析

在图2的全波列中，占主要地位的是弯曲波成分.为了分析弯曲波的频散特征，可以采用频域法从时域波形中提取各频率下弯曲波的慢度.本文利用频域加权相似法对图2中两个模型的全波波列进行处理，提取了偶极全波波列的弯曲波.

图3a中黑色和红色曲线分别表示模型1和模型2中第一道波列的频谱.图3b中，黑色和红色圆圈状曲线分别为对模型1和模型2的波列提取的弯曲波频散，对它们利用正切（或反正切）函数进行拟合，分别得到图中黑色和红色实线，利用频散方程［19］计算两个模型的弯曲波频散见三角状曲线，这三种方法求解的频散吻合性较好.拟合的好处是可以得到宽频带范围内的弯曲波，有效地消除了频散曲线中异常数据点对反演结果的影响，提高了反演计算的准确性.此外，从图3b中可以看出频率在1kHz以下时，红色和黑色实线几乎重合，并且随着频率向低频移动，二者的弯曲波慢度都趋向于模型1中地层横波慢度，这也就帮助确定发生径向变化前地层横波慢度的数值，在反演计算过程中，此数值为参考模型的地层横波慢度.

图2 偶极子全波列波形图，模型1（黑色）和模型2（红色）Fig.2 Dipole waveforms of model one（black）and model two（red）

3 反演方法

由图2和图3可以看出，地层横波慢度在径向上的变化分布影响了偶极子弯曲波传播特征，因此弯曲波传播特征的改变可以被用来确定地层横波慢度在径向上的变化情况，也即地层性质的改变.

本文选择的反演方法结合了微扰法和BG理论，根据此反演方法，待求模型的井外有一个横波慢度的径向变化带，此变化带是由于参考模型（参考模型的井外是无限大均匀地层）的横波慢度在近井壁处发生径向扰动所致.处理待求模型的偶极全波波列得到的弯曲波慢度在低频处的趋近值，可以认为是参考模型的地层横波慢度srefer.利用牛顿－拉夫森（Newton－Raphson）方法［20］计算频散方程：

可以求解出参考模型的弯曲波频散vrefer＝ω／k，进而求解出扰动部分贡献的弯曲波频散∂v＝vrealvrefer，其中vreal表示提取的弯曲波频散，符号∂及其后面符号代表小量的意思.利用∂v，得到地层横波慢度发生的扰动∂ss，它与参考模型的横波慢度之和，即为实际模型的横波慢度径向分布.

3.1 微扰法

在柱坐标系下，根据微扰法，当参考模型空间中有扰动发生后，扰动后的模型的弯曲波相速度也相应的发生了改变：

图3 （a）两个模型的第一道波列的频谱，（b）两种模型的弯曲波频散：提取（圆圈状）、拟合（线状）和利用频散方程计算（三角形）Fig.3 （a）The first spectrum of each model；（b）Flexural dispersion curves of each model：the extractive ones（circle），the fitted ones（line）and the calculated ones（triangle）based on dispersion equation

因此只要已知没有扰动时参考模型的相速度，就可求解发生扰动后的模型的相速度：

因此根据方程（9），当地层横波慢度从井壁r0处到径向半径b处，发生了扰动可以计算地层横波慢度从井壁r0处到b发生径向扰动后的弯曲波频散，只是这里的积分区域变成了r0到b，即

3.2 BG理论

Backus和Gilbert提出由于反问题解的非唯一性，不去求反问题解本身而去求它的某种平均乃是消除反问题非唯一性的一种途径，这种平均称为解估计.这里结合微扰法和BG理论，给出了利用弯曲波频散求解地层横波慢度的主要公式.

假设已知在频率点ωi处（i＝1，…，N），弯曲波相速度的相对变化∂vi／vi，是由r′处剪切模量发生的相对扰动∂μ（r′）／μ造成的，即有：

根据BG理论，在反演∂μ（r′）／μ时，不去求解∂μ（r′）／μ本身，而是求解∂μ（r′）／μ的估计∂μ＾（r′）／μ具有更好的稳定性.如下：

因此只要求解出待定系数组｛ai｝i＝1，…，N，就能计算出r′处的剪切模量的相对变化的估计.Backus和Gilbert的研究给出了系数组的一系列准则和方法，具体的过程见参考文献［5］.

下面介绍如何提高所采用的结合微扰法和BG理论的反演方法的计算效率，这也是本文的重点.

根据BG理论，正确评价解估计的准则是在展布函数和方差之间取合理的折衷.Backus和Gilbert建议采用以下目标函数：

其中，aTSa是平均核函数A（r，r′）的展布函数，aTEa是解估计的方差，β为折衷参数且β∈（0，π／2）.

求解方程（15）的过程是个对β不断寻优的过程，采用传统的迭代方式，计算效率较低，本文采用了1995年杨文采提出的改进措施［21］.在迭代开始时，先求低分辨率小方差的解估计，随非线性迭代逐步提高分辨率.具体地，在开始时取折衷因子β0＝π／2，随迭代逐步减小.即：

迭代次数q增加时，βq减小，这样的选取方式避免了内循环中解线性方程组的无谓计算，大大降低了计算成本.

3.3 反演流程

结合微扰法与BG理论，利用弯曲波频散反演地层横波慢度径向分布的基本流程可以归纳为：

（1）井孔半径、井孔流体参数和地层参数（除横波慢度外）已知，从波形数据中提取弯曲波频散vreal，对其进行拟合得到了宽频带范围的弯曲波慢度，其在低频处趋近的数值即为参考模型的横波慢度；

（2）对于选定的参考模型，计算偶极频散方程，得到参考模型的弯曲波频散vrefer；

（3）比较vreal与vrefer，若二者之间的差别小于某一小值min，则认为地层横波慢度没有发生径向扰动，否则认为地层横波慢度相对参考模型发生了径向扰动，则开始反演计算：

首先计算核函数（10）式，其次联合（10）—（15）式计算出待定系数组｛ai｝，i＝1，2，…，N，则得到了径向半径r′处的地层横波慢度的解估计：

改变r′值就可得到井外地层横波慢度的径向分布，它与地层横波慢度真实值之间相对误差为：

4 计算实例

本文假设模型2的偶极全波列为实际测得数据.为了分析结合微扰法和BG理论反演方法的特性，本文对偶极全波列存在以下5种情况时的数据进行了处理，得到了各情况下地层横波慢度的径向分布.这5种情况是：（1）无噪声和参数误差；（2）井孔流体慢度sf存在10%以内误差；（3）井外地层纵波慢度sp存在10%以内误差；（4）全波波列存在SNR＝20dB的噪声；（5）全波波列存在SNR＝10dB的噪声.

4.1 无噪声和参数误差

假设图2中的红色波列数据已知，其所对应的模型井外地层横波慢度径向分布情况未知.根据第2章介绍的反演方法，得到的参考模型的弯曲波频散见图4a中黑色曲线，处理图2中的红色波列得到的未知模型的弯曲波频散见图4a中红色曲线.本文选择6个频率点下两个模型的弯曲波慢度之差作为反演计算的输入参数，相邻点之间频率间隔为Δf＝1.0kHz，每个频率点下两个模型的弯曲波慢度用直线连接，见图中圆圈所标示.

图4b中红色曲线为反演计算得到的地层横波慢度的径向分布，黑色阶梯状曲线为地层横波慢度真实的径向分布.反演结果相对误差的径向分布见图4c，从图中可以看到反演结果相对误差基本控制在10%以内，并且离井壁越近，误差相对越大.

4.2 其他参数存在误差

实际数据比理论数据要更复杂，可能导致获取的其他地层参数或者井孔参数不准确，为此，我们还分析了这些参数存在误差时对反演结果的影响.本文首先研究了井孔流体慢度sf分别存在－10%、－5%、5%和10%的误差时对反演结果的影响.反演的与实际的地层横波慢度的径向分布见图5a，反演结果的相对误差的径向分布见图5b.我们注意到，sf存在误差与sf无误差的反演结果保持了很好的一致性.

其次研究了井外地层纵波慢度sp分别存在－10%、－5%、5%和10%的误差时对反演结果的影响.反演的与实际的地层横波慢度的径向分布见图6a，反演结果的相对误差的径向分布见图6b.我们注意到，sp存在误差与sp无误差的反演结果也保持了很好的一致性，这也是计算时假设地层纵波慢度为一固定值的原因.

4.3 存在噪声

实际数据中往往存在噪声，为此还应针对波列数据存在噪声的情况做进一步的分析.这部分的工作是通过对图2中的红色波列分别添加SNR＝20dB和SNR＝10dB的高斯白噪声各4次来展开的.

图7显示了分别添加了一次SNR＝20dB和SNR＝10dB的高斯白噪声后的波形曲线和频谱曲线.图7a中黑色、红色和蓝色曲线分别表示：无噪声时、存在SNR＝20dB和SNR＝10dB的噪声时的偶极全波列波形；图7b的曲线为相应的第一道波形，图7c为相应的第一道波形的频谱（图中标示了果的准确性大大降低.拟合后的弯曲波频带范围变宽，可利用的弯曲波数据增加，进而降低了噪声的干扰.加入4次SNR＝20dB噪声后反演的与实际的地层横波慢度的径向分布见图8e，反演结果相对误差的径向分布见图8f.加入4次SNR＝10dB噪声后反演的与实际的地层横波慢度的径向分布见图9e，反演结果相对误差的径向分布见图9f.从图8f和9f可见：存在SNR＝20dB和SNR＝10dB噪声时，相对误差也基本控制在10%以内，并且与没有噪声情况下的误差吻合度较高，这说明采用的反演方法对噪声不敏感.

图7 （a）偶极全波列波形；（b）第一道波列；（c）第一道波列的频谱；（d）SNR＝20dB时八道波列的频谱；（e）SNR＝10dB时八道波列的频谱Fig.7 （a）Dipole full waveforms；（b）The first waveform of each case；（c）The first spectrum of each case；（d）Eight spectrums when SNR＝20dB；（e）Eight spectrums when SNR＝10dB

5 结论与讨论

本文以井外地层横波慢度径向分层模型为例，计算了模型的偶极全波波列，提取了模型的弯曲波频散，利用此弯曲波慢度与井外地层均匀时弯曲波慢度在多频率点下的差异，采用结合微扰法和BG理论的反演方法，得到了径向分层模型井外地层横波慢度的径向分布，并且讨论了如何提高反演计算的效率，以适应现场对数据处理速度的要求.

此外，为了分析采用的反演方法的特性，本文分别对以下5种情况的偶极全波波列数据进行处理，得到了各情况下地层横波慢度的径向分布.这5种情况是：（1）无噪声和参数误差时；（2）井孔流体慢度sf存在10%以内误差时；（3）井外地层纵波慢度sp存在10%以内误差时；（4）全波波列存在SNR＝20dB的噪声时；（5）全波波列存在SNR＝10dB的噪声时.结果表明：（2）和（3）的反演结果与（1）的反演结果保持了很好的一致性，这是因为sf和sp对弯曲波的影响较小的缘故；在（4）和（5）两种情况下，通过对提取的弯曲波频散进行拟合，得到了较宽频带范围内的弯曲波数据，有效地消除了频散曲线中异常数据点对反演结果的影响，使得（4）和（5）的反演结果相对于（1）的反演结果都没有发生明显的改变，这说明采用的反演方法对噪声不敏感.

图9 （a—d）SNR＝10dB时提取（圈状）和拟合（线状）的弯曲波频散；（e）反演的与实际的地层横波慢度的径向分布；（f）反演的与实际的地层横波慢度之间相对误差的径向分布Fig.9 （a—d）Extractive（circle）and fitted（line）flexural dispersion curves when SNR＝20dB；（e）Inverted and actual radial profiles of shear wave slowness；（f）Radial profile of relative error between inverted and actual shear wave slowness

以上反演结果都证实了，结合微扰法和BG理论的反演方法具有很好的鲁棒性，可以被用于反演井壁周围地层横波慢度的径向分布，进而了解井壁周围地层的性质.

（

）

［1］ Schmitt D P.Shear wave logging in elastic formations.JournaloftheAcousticalSocietyofAmerica，1988，84（6）：2215－2229.

［2］ Winbow G A.A theoretical study of acoustic S－wave and Pwave velocity logging with conventional and dipole sources in soft formations.Geophysics，1988，53（10）：1334－1342.

［3］ Ellefsen K J，Cheng C H，Toksöz M N.Applications of perturbation theory to acoustic logging.JournalofGeophysical Research，1991，96（B1）：537－549.

［4］ Backus G，Gilbert F.Uniqueness in the inversion of inaccurate gross earth data.Philosophical Transactions of the Royal Society A：Mathematical.PhysicalandEngineering Sciences，1970，266（1173）：123－192.

［5］ Burridge R，Sinha B K.Inversion for formation shear modulus and radial depth of investigation using borehole flexural waves.SEG Expanded Abstracts，1996，158－161.

［6］ Sinha B K，Burridge R.Radial profiling of formation shear velocity from borehole flexural dispersions.IEEE UltrasonicsSymposiumProceedings，2001，391－396.

［7］ Sinha B K，Asvadurov S.Dispersion and radial depth of investigation of borehole modes.GeophysicalProspecting，2004，52（4）：271－286.

［8］ Sinha B K，Vissapragada B，Renlie L，et al.Radial profiling of the three formation shear moduli and its application to well completions.Geophysics，2006，71（6）：E65－E77.

［9］ Braunisch H，Habashy T M，Sinha B K，et al.Inversion of borehole dispersions for formation elastic moduli.IEEE UltrasonicsSymposium，2000，551－556.

［10］ Tang X M，Patterson D J.Mapping formation radial shearwave velocity variation by a constrained inversion of borehole flexural－wave dispersion data.Geophysics，2010，75（6）：E183－E190.

［11］ Yang J，Sinha B K，Habashy T M.A Parameterized－modelbased radial profiling for formation shear slowness in cased boreholes.SEG San Antonio 2011Annual Meeting，2011，449－453.

［12］ Yang J Q，Sinha B K，Habashy T M.A theoretical study on formation shear radial profiling in well－bonded cased boreholes using sonic dispersion data based on a parameterized perturbation model.Geophysics，2012，77（3）：WA 197－WA 210.

［13］ 陈雪莲，王瑞甲.横向各向同性弹性地层井孔中模式波的探测深度.吉林大学学报（地球科学版），2008，38（3）：502－507.

Chen X L，Wang R J.Investigating depth of mode waves in the borehole surrounded by transversely isotropic elastic formation.JournalofJilinUniversity（EarthScience Edition）（in Chinese），2008，38（3）：502－507.

［14］ Nolte B，Huang X J.Dispersion analysis of split flexural waves：Borehole acoustics and logging and reservoir delineation consortia annual report.Massachusetts Institute of Technology，1997.

［15］ Tang X M，Li C，Patterson D J.A curve－fitting technique for determining dispersion characteristics of guided elastic waves.Geophysics，2010，75（3）：E153－E160.

［16］ Paillet F L，Cheng C H.Acoustic waves in boreholes.Boca Raton／Ann Arbor／Boston London：CRC Press，1991.

［17］ Schmitt D P.Effects of radial layering when logging in saturated porous formations.JournaloftheAcoustical SocietyofAmerica，1988，84（6）：2200－2214.

［18］ 张海澜，王秀明，张碧星.井孔的声场和波.北京：科学出版社，2004.

Zhang H L，Wang X M，Zhang B X.Acoustic wavefield and waves in the borehole（in Chinese）.Beijing：Science Press，2004.

［19］ Tang X，Wang T，Patterson D.Multipole acoustic loggingwhile－drilling.／／72nd Annual International Meeting，SEG Expanded Abstracts，2002，364－368.

［20］ Press W H，Teukolsky S A，Vetterling W T，et al.Numerical recipes in Fortran 77：The art of scientific computing.Cambridge：Cambridge University Press，1992.

［21］ 杨文采.地震波场反演的BG一逆散射方法.地球物理学报，1995，38（3）：358－366.

Yang W C.BG－inverse scattering method for inversion of seismic wavefield.ChineseJ.Geophys.（in Chinese），1995，38（3）：358－366.