0 引 言

在多波多分量地震勘探中，当转换横波斜交穿过由裂缝引起的各向异性介质时会发生分裂，产生偏振方向平行裂缝走向的快横波和垂直裂缝走向的慢横波，我们称之为横波分裂现象.分裂的快慢横波的强弱与裂缝的强度密切相关(缪林昌，1994).通过对横波分裂现象的分析可以得到裂缝的发育方向和密度等参数，利用这些参数可以更为准确地预测裂缝.其中，快横波的偏振方位与裂缝方位走向密切相关，因此我们可以将求取裂缝走向的问题转化为求取快横波的偏振方位角度的问题；又由于快、慢横波的传播时差反映了裂缝发育密度，所以可以将裂缝发育密度的预测问题转化为求取快、慢横波的传播时差问题.利用多分量地震中转换横波分裂这一特性研究裂缝方向及其发育程度是目前最直接、最可靠裂缝检测方法(李彦鹏和马在田，2000；张晓斌等，2003；郭桂红等，2008；徐天吉等，2011).

粒子群算法(Particle Swarm Optimization，PSO)是20世纪90年代发展起来的一种新型仿生模拟进化算法.该算法的思想就是在可行解空间中随机初始化一群粒子，每个粒子都为优化问题的一个可行解，并由目标函数为之确定一个适应度值，同时还有一个速度来决定其方向和距离.通常粒子将追随当前的最优粒子而动，并经逐代搜索最后得到最优解(Kennedy and Eberhart，1995).由于粒子群算法本身简单易于实现、参数设置少，可以解决复杂的非线性问题，在短期内得到很大发展，己在函数寻优、神经网络训练、数字电路优化、问题以及工程应用等多个优化及反演领域得到了成功的应用(朱童等，2011).模拟退火法(Simulated Annealing Method，SA)源于统计热力物理学，Dariu于2005年利用模拟退火原理提出了一种全局最优化算法来估计横波分裂属性，研制了一种自动逐层估算极化方向和快慢波时间延迟技术(Dariu et al.，2005).模拟退火法的基本思想是生成一系列参数向量模拟粒子的热运动，通过缓慢地减小一个模拟温度的控制参数，使模拟的热系统最终冷却结晶达到系统能量最小值(Kirkpatrick et al.，1983).

在本文中，针对粒子群优化算法和模拟退火算法各自的优势与不足，为了增强粒子群优化算法的全局搜索能力和跳出局部最优的能力，增加粒子群的多样性，提高算法的收敛速度，笔者将模拟退火法与收缩因子粒子群算法进行融合，结合二者的优点，利用快慢横波的Pearson相关系数函数作为目标函数，来自动求取裂缝的方位与密度，实现快慢横波的分离，达到了裂缝属性识别的效果.

1 算法原理 1.1 快慢横波Pearson相关系数

转换横波在由裂缝引起的各向异性介质中传播时，会产生横波分裂，转换横波被分裂为快横波S₁(t)和慢横波S₂(t).在观测坐标系统上，快慢横波在x方向的投影和构成径向分量S_x(t)，在y方向的投影和构成横向分量S_y(t)，即(Alford，1986)

(1)

经坐标变换可求得快、慢波分量S₁(t)和S₂(t)的公式为

(2)

其中，θ为测线x方向与裂缝走向的夹角.

快、慢横波是经由同一个横波分裂而得出来的，只是存在时差上的差异而已，因此二者的波形应该具有最大的相似程度.而数学中的Pearson相关系数主要用来描述两个向量上的数据是否在一条线上，我们可以用其来对快慢波波形的相似性进行判别.其数学表达式为(张宇镭等，2005)：

(3)

其中，N为向量元素个数，X、Y为两个等长向量.

首先选取采样点为N的两个地震记录S_x(t)和S_y(t)，然后给定一个角度θ，利用公式(2)可得到快、慢横波S₁(t)和S₂(t)，再给定一个时差Δτ，令X=S₁(t+Δτ)，Y=S₂(t)，将X、Y代入公式(3)就得到快、慢横波的Pearson\相关系数公式为(王凯等，2012)

(4)

1.2 粒子群算法

粒子群算法(PSO)由Kennedy和Eberhart共同提出的，其基本思想是通过对鸟类的群体行为进行建模与仿真研究结果而得到启发(Kennedy and Eberhart，1995).此后，Shi Y和Eberhart又提出了惯性权重这一概念，并将其引入到基本粒子群算法的速度进化公式中(Shi and Eberhart，1998，2001；Eberhart and Kennedy，2001)，并提出了在进化过程中动态的调整惯性权重来平衡收敛的速度，该进化方程被称为标准粒子群算法(Shi and Eberhart，1998).2002年，Clerc等在标准粒子群算法的基础上提出的带收缩因子的粒子群算法，与标准粒子群优化算法相比，收缩因子粒子群优化算法有更快的收敛速度，跳出局部最优解的能力也有所提高(Clerc and Kennedy，2002).其位置和速度更新方程为

(5)

(6)

上式中，下标i表示第i个粒子；下标m表示粒子的第m维；常数c₁和c₂为学习因子；χ=2/2-φ-(φ²-4φ)^1/2称为收缩因子，且φ=c₁+c₂，φ>4；r₁和r₂是取值范围为[0，1]之间的随机数；v_im^k和v_im^k+1分别表示第k次和第k+1次迭代时第i个粒子的速度矢量的第m维分量；x_im^k和x_im^k+1分别表示第次和第次迭代时第i个粒子的位置矢量的第m维分量；pbest_im是第i个粒子的个体最佳位置pbest的第m维分量gbest_m是种群最好位置gbest的第m维分量.

运用粒子群算法进行裂缝属性识别就是将裂缝方位角θ和快慢波时差Δτ作为一个二维粒子，即x_i=(θ_i，Δτ_i)，把快慢横波的Pearson相关系数函数作为适应度函数，求解Pearson相关系数最优解的过程.

1.3 模拟退火粒子群算法

模拟退火粒子群算法(SAPSO)以粒子群算法为主体流程，这样保证了新算法继承了PSO算法简单，容易实现的特性.然后引进模拟退火，在每个粒子的速度和位置更新过程中加入Metropolis准则，对PSO进化后的适应度值按Metropolis准则接受优化解的同时以一定的概率接受恶化解，算法从局部极值区域中跳出，然后调整温度使其逐渐下降，粒子便逐渐形成低能量基态，收敛至全局最优解.

在引入模拟退火的过程中，温度更新函数为

(7)

其中，k是一个接近1的常数，一般取0.5~0.99，该衰减函数对控制参数的衰减量是随算法迭代次数递减的.

基于模拟退火粒子群算法的的裂缝属性识别方法的执行流程如下：

(1) 首先在径向分量x和横向分量y地震记录上，分别选取采样点为N的两个记录S_x(t)和S_y(t)作为裂缝属性识别的数据.

(2) 随机选取N个二维粒子(θ，Δτ)，初始化每个粒子的位置和速度，则第i个粒子的位置和速度为x_i=(θ_i，Δτ_i)，v_i=(v₁，v₂).设定最大迭代次数k_max和学习因子c₁和c₂，设置初始温度T₀，其中T₀应该充分大，以及允许变坏的范围Fval

(3) 运用公式(4)计算每个粒子对应的适应度.在公式(5)中pbesti是第i个粒子的个体最佳位置； gbest是所有粒子的全局最优位置.

若f(pbest_i)≤f(x_i)，则更新pbest_i的位置，即pbest_i=x_i；

若f(gbest)≤f(pbest_i)，则更新gbest的位置，即gbest=pbest_i.

(4) 根据公式(5)和(6)更新每个粒子的位置和速度.

考虑位置：若x_i^k+1>x_max，则x_i^k+1=x_max;若x_i^k+1<-x_max，则x_i^k+1=-x_max;否则x_i^k+1不变.

考虑速度：若v_i^k+1>v_max，则v_i^k+1=v_max;若v_i^k+1<-v_max，则v_i^k+1=-v_max;否则v_i^k+1不变.其中，x_max、v_max分别为最大位置和最大速度，都是常数.

(5) 计算第i粒子两次位置变化适应度变化量，Δf=f(x_i^k+1)-f(x^k_i)，如果Δf≤Fval，则接受粒子的新速度v_i^k+1为当前速度，接受新产生的位置x_i^k+1为当前最优点；否则，计算p=exp(-Δf(T)，如果p>rand(0，1)，则接受粒子的新速度v_i^k+1为当前速度，接受新产生的位置x_i^k+1为当前最优点；否则，v_i^k、x_i^k不改变.

(6) 如果迭代次数达到最大迭代次数k_max，算法收敛，最后一次迭代的全局最优值gbest就是我们所求的最优值点，所对应的粒子位置x_max=(θ_max，Δτ_max)就是获得的最优解；否则，根据公式(7)重新更新退火速度T，更新迭代次数k=k+1，返回步骤(3)，算法继续运算.

(7) 步骤(6)中求取的θ_max和Δτ_max即为裂缝方位角和快慢波时差，最后根据Thomsen参数与时差Δτ_max和裂缝密度e的关系式，即可求解裂缝密度e.

2 模型试算

在本文中，模型试算数据选取的是二维三分量的EDA介质正演模拟数据.介质模型大小为400×400个网格点，Δt=1 ms，Δx=Δy=Δz=5 m，子波类型为雷克子波，主频为30 Hz，纵波速度为2800 m/s，横波速度为1800 m/s，模型的裂缝密度为0.1，方位角为45°.裂缝属性识别方法就是利用模拟退火粒子群算法自动求取裂缝角度和密度，其中取c₁=c₂=2.05，粒子数目为10，迭代次数为30次.初始温度T₀=f(gbest)/ln5.

2.1 单道地震记录的裂缝属性识别

在正演模拟数据中随机选取一道地震记录作为实验数据.根据模拟退火粒子群算法求取裂缝属性的步骤分别进行了单次运行计算和多次运行计算.单次运行计算中，得到的方位角为47°，误差为4.444%；时差为17，裂缝密度为0.1024，密度误差为2.395%.图 1中显示的是显示的是快、慢横波Pearson相关系数的迭代结果，可以看出SA-SFPSO算法的收敛速度很快.图 2显示的是通过求取出来的裂缝方位角和时差，对快慢横波进行了分离.为了验证方法的稳定性，又进行了100次的运行，并将得到的结果与方差与基本粒子群算法相比较，如表 1所示.

此外，实验为验证算法跳出局部最优解的能力，还进行了求1000次单独运行中出现的适应度低于0.9的次数，为准确考虑，对实验进行100次独立运算，求取其结果均值和方差，并将结果与基本粒子群算法进行了比较，如表 2所示.

2.2 多道地震记录的裂缝属性识别

从地震正演模拟记录中选取100道地震记录作为实验数据，为了验证算法的抗噪性能，同时加入10 dB随机高斯白噪声.图 3显示的是信噪比为10 dB的地震记录，图 4表示的是裂缝方位角度玫瑰图，图中轴的长度表示识别出的对应角度的数目，轴的长度越长表示识别出该角度的数量越多，从图 4中我们可以看出长轴大致分布在45°左右，而单次运算的角度均值为46.07°.图 6表示的各道的快慢波时差.图 5显示的是裂缝密度直方图，该图中轴长同样表示对应裂缝密度的数目，轴长越长表示识别出对应裂缝密度的数量越多，图中显示长轴均分布在0.1附近，单次运算的裂缝密度均值为0.106.同样为了验证算法在多道数据处理中的稳定性，进行了100次运行，并得到的结果与方差与基本粒子群算法的相比较，如表 3所示.

3 结 论 3.1  

本文中，为了克服传统粒子群算法的缺点，提高算法的性能，笔者首次将模拟退火法与收缩因子粒子群算法进行融合，将模拟退火算法的概率突跳能力引入到收缩因子粒子群优化算法，让新算法在具有较快收敛速度的同时还具有较强的全局寻优能力，并首次利用该融合算法与Pearson相关系数法相结合自动识别裂缝方位角及密度，模型试算的结果表明：

(1)在单道地震数据处理中，本方法能够可以很快的收敛于全局最优解，而且运算结果的稳定性与精度都非常高，快慢波也能得到很好的分离.

(2) 为了验证在噪声环境中本方法的有效性，在处理多道实验数据时加入了随机噪声，从运算结果分析比较可得，本方法得到的结果也相当精确，抗噪性能较强.

(3) 基于模拟退火粒子群算法与Pearson相关系数法相结合的方法是一个有效的识别裂缝属性的方法.该算法具有两种算法的优点，有良好的跳出局部最优解的能力，能够更快的收敛于全局最优解，而且稳定性与精度都较高.

3.2  

虽然本文提出的方法在合成地震数据处理中能取得很好的应用效果，并在一定程度上提高了裂缝属性识别的稳定性和精确性，但是由于缺乏实际地震数据处理，笔者并未验证该算法在实际数据中的应用效果，所以下一步工作就是将该算法应用到裂缝性油气藏等实际的地震数据处理中.而且本文算法目前只适应于单一裂缝地层，但实际地层一般会包含很多裂缝地层，因此可以将本算法与层剥离技术结合起来求取各层的裂缝属性.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持！