1 引言

对套管井外地层进行有效的评价在油气藏开发的中后期非常重要，套管外地层受异常地应力的影响而遭到损坏(李刚，2005; 唐晓明等，2015；陈雪莲等，2015；曹景记等，2014)，导致地层在径向上表现出非均质性，在声波测井中即为地层速度的径向变化.确定套管井外地层速度径向变化剖面对井眼稳定性评价及获取有效的生产动态监测数据具有重要的意义.近年来，国内外学者对利用声波测井数据反演裸眼井地层横波速度径向分布进行了大量的研究.Sinha等提出了结合微扰法和Backus-Gilbert理论的反演算法(Bruridge and Sinha, 1966；Backus and Gilbert, 1970；Sinha and Bruridge, 2011；Sinha and Ascadurov, 2004；马明明等，2013)，微扰法精度较低，而且逐点反演地层横波速度的径向分布，计算效率也较低.Tang和Patterson (2010)建立了参量化的、高频约束条件下的反演方法，在反演的过程中对弯曲波进行高频约束，但高频弯曲波的激发强度较小，易受噪声的干扰.Yang等(2012)建立了套管井模型的参量化的目标函数，假定地层剪切模量是径向距离r的指数函数，将横波速度径向分布反演问题转换成三个控制参数的反演，用高斯牛顿法求解了关于这三个参数的非线性方程组.

本文基于Yang等(2012)建立的参量化的目标函数，选择双层参考模型，采用修正的微扰法计算了套管井地层横波速度径向分层模型的弯曲波频散曲线，验证了该计算方法的准确性.采用高斯牛顿(Gill et al., 1981)和快速模拟退火法(何峰江和陶果，2005)求解反演目标函数，得到了地层剪切模量参数，进而得到地层横波速度的径向分布，计算结果与理论模型误差较小.分析了偶极弯曲波频段、套管横波速度对反演结果的影响，发现采用偶极弯曲波激发强度较高的频段与采用全频段的反演结果相近；套管的横波速度准确度越高，反演结果越准确.对采用高斯牛顿法和快速模拟退火法的计算结果进行对比发现：两种方法有相近的计算精度；高斯牛顿法具有较高的计算效率，对反演参数初始值的选取要求较高；快速模拟退火方法的计算效率略低于高斯牛顿方法，但受初始值的选取影响较小，实际处理中应选择快速模拟退火法.针对地层横波速度径向分层模型，理论模拟了不同中心频率声源的偶极全波列，对理论模拟数据进行处理得到了地层横波速度的径向分布，结果与模型一致.此外，通过对结果分析还得出反演结果受声源中心频率影响较小.

2 套管井中偶极模式波的频散方程

建立同心圆柱层状模型模拟套管井几何结构，如图 1所示.从内到外依次为井内流体、套管、速度变化带、原状地层.此模型下，研究的问题为同心圆柱层介质在轴向中的多极声场分布和声波传播问题.井内、外的位移势在波数-频率域(时间因子e^iωt)的表达式为(Schmitt，1988)：

(1)

其中a表示套管内径，K_n表示第二类贝塞尔函数，代表由内向外传播的声波，I_n表示第一类变型贝塞尔函数，表示由外向内传播的声波. , 和分别是流体纵波、固体纵波和固体横波径向波数，v_f、v_p和v_s分别是流体纵波、固体纵波和固体横波速度，k_f=ω/v_f，k_p=ω/v_p和k_s=ω/_s分别是对应波数.n=1代表偶极声源，式中A′_n、A_n、B_n、C_n、D_n、E_n和F_n是待定系数，由井壁及套管内外的边界条件确定.柱坐标系下位移分量的表达式为：

(2)

利用位移分量，可得柱坐标系下的应变分量：

(3)

由应变分量，利用胡克定律，可得到应力分量：

(4)

其中λ为拉梅系数、μ为剪切模量，对于井内流体，取μ=0，e=(e_rr+e_θθ+e_zz)是体应变；i=j时δ_ij=1否则δ_ij=0，下标i，j分别代表r，θ，z，在r=a处，径向位移和径向正应力连续，切向位移为零，如(5)式：

(5)

如果在半径r_m(下标m表示径向的第m层)处的界面两侧都是固体，在界面上胶结良好，则界面上的位移和应力分量连续，即

(6)

层与层之间的声场由汤姆森-哈斯克(Thomson-Hanskell)传播矩阵来连接，反复利用各层介质的传递矩阵和边界条件，得到如下的矩阵方程：

(7)

解出A′_n、B_n、D_n和F_n，利用递推关系可以求出全部的系数，由克莱姆法则

(8)

式中det表示矩阵的行列式，矩阵M^ij是余矩阵，是矩阵M去掉第i行、第j列得到的，式(8)中的极点由(9)式确定

(9)

方程(9)称为频散方程.该方程是相对于k的非线性方程，常用牛顿-拉夫森方法来进行求解，得到弯曲波相速度随频率变化的方程.

(10)

3 套管井中偶极模式波频散分析

套管井高速地层的偶极弯曲波已被深入研究过，套管的存在对地层参数反演的影响研究较少，尤其是套管外存在地层变化带的情况.本小节分析了套管参数对弯曲波频散的影响，变化程度用灵敏度sen_q表示，定义如下：

(11)

下标q表示某参数.

裸眼井半径为0.080 m，地层参数见表 1.套管井模型如图 1所示，套管与地层胶结良好，套管内径a=0.080 m，套管外径b=0.088 m，变化带外半径c=0.3 m，变化带参数除地层横波速度外其他参数与原状地层相同，具体见表 1，变化带横波速度随半径r的变化见表 2.

图 2a为裸眼井和套管井的偶极弯曲波频散曲线的理论计算结果.图 2b是裸眼井和套管井的偶极弯曲波的激发强度曲线.图 2c是偶极弯曲波频散对套管横波速度、套管纵波速度、套管厚度(变化套管外半径)、套管密度的灵敏度曲线.

由图 2c可以看出，套管纵波速度和密度对偶极弯曲波频散的影响非常小，套管厚度对偶极弯曲波频散的影响稍大，套管横波速度对偶极弯曲波频散的影响是最大的.因此利用弯曲波频散反演套管井地层横波速度径向分布时，必须考虑套管横波速度的影响.

4 套管井外地层横波速度径向分布反演

对于套管井外地层横波速度径向分层模型，套管参数已知时，可根据式(9)来建立反演目标函数：

(12)

其中N表示选择的频率点数，m为要反演的层数，要反演每个层的厚度和速度两个参数，至少要保证N≥2m.实际计算时，有2m个参数需要确定，由于反演计算的多解性，直接求解(12)式是不可行的.井外地层变化带的剪切模量可以表示为(Yang et al., 2012)：

(13)

其中，c_μ是变化带外径，μ_F是原状地层剪切模量，Δμ_A是变化带剪切模量相对于原状地层剪切模量的最大扰动，σ(σ>0)是指数函数的衰减系数.根据(13)式，已知变化带参考剪切模量为μ_A时，相应的扰动可表示为

(14)

参数c_μ、Δμ_A和σ的不同组合即可以用来表征不同的横波速度径向分布.有效的反演算法基于准确的正演模拟结果，接下来对套管井外地层横波速度径向分层模型的频散方法进行分析，目前常用的为微扰法.

4.1 微扰法计算偶极弯曲波频散 4.1.1 传统计算方式

在柱坐标系下，根据微扰法的计算方式，当模型空间中有扰动发生后，模型的相速度也相应的发生变化：

(15)

其中c_ijkl=λδ_ijδ_kl+μ(δ_ikδ_jl+δ_ilδ_jk)为弹性系数，，λ和μ为拉梅系数，*表示复共轭，只有剪切模量发生扰动时，根据应变与位移的关系，有：

(16)

其中，

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

(17)

4.1.2 修正计算方式

一种修正的计算方法是：定义变化带剪切模量μ_A发生的平均扰动为Δμ_A，根据灵敏度的定义式(11)有：

将(15)和(16)式代入该式得到：

得到等效的变化带横波速度

(18)

可得变化带地层剪切模量发生的平均扰动为：

等效的变化带地层横波速度是：

(19)

采用等效横波速度重新建立井孔模型，此时v_sA和v_sF是波数k的函数，利用传播矩阵，计算频散方程(9)可以得到扰动发生后的角频率ω_new，根据得到新的弯曲波频散曲线.

(20)

4.1.3 偶极弯曲波频散计算结果分析

由图 2c可知，套管横波速度对偶极弯曲波频散影响较大，因此利用弯曲波频散反演地层横波速度径向分布时，必须考虑套管横波速度的影响.本文选择参考模型时假设套管的参数已知，将参考模型分单层和双层两种情况.单层参考模型的套管外有一个无限大均匀弹性地层，与套管胶结良好，井孔和套管参数见表 1，地层横波速度2500 m·s^-1.双层参考模型的套管外有一个变化带，与套管和原状地层胶结良好，变化带横波速度为2000 m·s^-1，单、双层参考模型和地层模型的横波速度径向分布见图 3a.图 3b为利用微扰法和传播矩阵法计算的弯曲波频散，图 3c为两种方法计算结果的相对误差.图中线1和线2分别表示选择单层参考模型时，微扰法传统计算方式和改进方式计算的结果，线3和线4分别表示选择双层参考模型时，微扰法传统计算方式和改进方式计算的结果，线5表示利用传播矩阵计算的弯曲波频散，为准确值.由图 3c可知采用改进的微扰法的计算结果更接近真实值，双层参考模型的计算结果优于单层参考模型，此外所有方法的相对误差受频率的影响比较明显，在高频时相对误差都趋近于固定值.

4.2 地层横波速度径向分布反演

由前文分析结果可知，采用双层参考模型及修正的微扰法计算的弯曲波频散与理论计算结果最为接近；不同频率下，地层横波速度对弯曲波频散的影响不同，因此可选择多频率点反演地层横波速度的径向分布.给定M组频散曲线的组合(ω_i, k_zi)，i=1, 2, …, M，根据方程(20)建立目标函数如下：

(21)

地层横波速度径向分布的反演问题，转化成为求解c_μ、Δμ_A和σ的非线性最小二乘问题.

4.2.1 反演计算流程

当M≥3时，可同时求解三个未知参数c_μ、Δμ_A和σ，具体的计算过程为：

首先，选择频率段为fre=(f_start, f_end)，频率采样间隔为Δfre，Δω=2πΔfre，得到组的弯曲波频散的组合(ω_i, k_zi)，i=1, 2, …, M；

其次，已知原状地层剪切模量μ_F，将地层剪切模量的径向分布(13)式和径向变化(14)式代入方程(18)，得到变化带的等效剪切模量(横波速度)Δμ_A；

最后，将(ω_i, k_zi)和Δμ_A代入(21)式后，利用反演方法求解c_μ、Δμ_A和σ.

4.2.2 反演结果影响因素分析

(1)频率范围

为了考察频率范围对反演结果的影响，首先选取有效激发的整个频率段的数据进行计算.由图 2b可知该频率范围为fre=(2.5 kHz, 13.5 kHz).分别采用高斯牛顿法和快速模拟退火法进行反演，计算结果如图 4所示.其中图 4a中点线和虚线分别是利用高斯牛顿法和快速模拟退火法得到的地层横波速度的径向分布；图 4b虚线和点线分别为两种方法的相对误差，由图可见两种方法的计算误差都很小，小于4%，说明选择整个频率段数据进行反演计算时，两种方法都是有效的，且具有基本相同的计算精度.

由图 2b可知，本文所采用套管井模型的偶极弯曲波的激发强度，在频段(2 kHz, 6 kHz)和(11 kHz, 14 kHz)范围内较小，容易受其他因素影响.频段(6 kHz, 11 kHz)内弯曲波的激发强度较大，故采用频段范围fre=(6 kHz, 11 kHz)进行计算，计算结果如图 5所示.其中图 5a中的点线和虚线分别为利用高斯牛顿法和快速模拟退火法得到的地层横波速度的径向分布；图 5b虚线和点线为两种方法的相对误差.由图可见两种方法的计算结果的误差都很小，都小于4%，说明选择该频段的数据即可进行有效的反演.

(2)反演参数初始值

求解非线性最小二乘方程组(21)的过程中，发现高斯牛顿方法的收敛性对c_μ、Δμ_A和σ的初始值依赖较高，具体见表 3，三个参数中任何一个参数的初始值发生10%以上的扰动时，都会导致迭代不收敛.采用快速模拟退火方法在计算成本上略有增加，但其收敛性受初始值的影响较小，具体见表 4.对比表 3和表 4可知，快速模拟退火法对初始值的依赖较小，更适用于实际数据的反演处理.

(3)套管横波速度

本文计算了套管横波速度误差分别为-10%、-5%、5%、0和10%时的反演结果.图 6a和图 7a分别为采用高斯牛顿法和快速模拟退火法的计算结果.图 6b和图 7b分别为反演结果与真实地层横波速度径向分布的相对误差.从相对误差分布图上可见，套管横波速度误差越大，反演结果的误差在近井壁处越大，并且当套管横波速度的误差为正时(比真实值大)，近井壁处的反演误差向负方向增大，当套管横波速度的误差为负时(比真实值小)，近井壁处的反演误差向正方向增大.

5 结论

本文采用双层参考模型，应用修正的微扰法得到了准确的弯曲波频散曲线.分别采用高斯牛顿法与快速模拟退火进行反演目标函数求解，得到了横波速度径向剖面的分布.通过模型试算，反演结果与理论模型一致.经过对比分析得到如下结论：

(1)利用激发强度较高的频段内的弯曲波数据进行反演，与利用全频段范围的弯曲波数据进行反演无明显的差异，在实际应用中可选择激发强度较高的频段内的数据进行处理，以提高信噪比.

(2)套管参数对反演结果的影响不可忽略，套管横波速度的影响最大，套管横波速度误差越大，近井壁处的地层横波速度反演结果误差越大.

(3)高斯牛顿方法的收敛性对初始值依赖性较大，而快速模拟退火方法对初始值依赖性较小，在实际应用中应选择快速模拟退火法进行反演.