0 引言

多分量叠前联合反演是一种综合利用纵波、横波以及转换波等地震记录中所含地下地层反射信息来预测介质弹性参数和物性参数(纵横波速度、密度和各向异性参数等)的有效地球物理方法(Ostrander，1984；张广智等，2014).目前，该方法在储层评价和油气预测等方面有着广泛的应用.随着油气勘探开发进程的逐步发展，各向异性理论研究以及复杂岩性油气藏勘探逐渐得到了人们广泛的关注，因此多分量叠前联合反演方法也由各向同性介质的三参数反演逐渐向各向异性介质的多参数反演方向发展.VTI介质就是一种常见的各向异性介质(如页岩等)，它是一种具有垂向同相轴的横向各向同性介质，是目前各向异性介质研究的热点(李勤等，2014).

各向异性的研究最早始于19世纪50年代，但是直到20世纪80年代才得到空前的发展，期间针对不同类型的各向异性介质提出了许多目前仍在广泛应用的理论.在这其中，横向各向同性介质(TI介质)最受人们的关注，其具体又可分为水平横向各向同性介质(HTI)和垂直横向各向同性介质(VTI).因为这两种介质都具有一个旋转对称轴，因此均可用五个独立的弹性参数和密度进行描述，这五个弹性参数被称为弹性刚度系数，可以用C₁₁，C₁₃，C₃₃，C₄₄和C₆₆进行表示.但是与VTI介质不同的是在HTI介质中弹性波的传播特征与方位角有关，为了区分这一差别，Thomsen等学者又将HTI介质称为方位各向异性介质或EDA介质，而用横向各向同性介质来专指VTI介质.同时，Thomsen(1986)针对横向各向同性介质提出了垂向纵横波速度(V_P0，V_S0)、密度(ρ)和三个弱各向异性参数(ε，δ和γ)来描述弹性波的传播特征.

在VTI介质的正演模拟方面，Daley和Hron(1977)根据界面两侧位移和应力连续条件，最先提出了利用相角表示的横向各向同性介质的反射和透射系数精确公式；Banik(1987)进一步解释了Thomsen各向异性参数的物理意义，并推导了VTI介质的PP波反射系数近似公式；Graebner(1992)基于特征值理论推导了用水平慢度表示的VTI介质的反射和透射系数精确公式；Rüger(1997, 1998)推导了VTI介质的PP波和PSV波反射系数近似公式.与Zoeppritz方程及其近似公式类似，上述这些公式仅仅反映了反射系数随界面两侧弹性参数和各向异性参数的变化，未考虑几何扩散、吸收衰减、透射损失和多次波等地震波传播效应.因此，许多学者(Kennett，1983；Fryer and Frazer, 1984；Mallick and Frazer, 1987；寻浩等，1997)先后研究了各向异性介质的反射率法，该方法能够精确地模拟水平层状介质的全波场响应.

在多分量叠前联合反演方面，前人也做了大量的研究.Stewart(1990)最早提出了实际的纵、横波联合反演方法，随后国内外许多学者针对各向同性介质对其进行了大量的研究(Fatti et al., 1994；Veire and Landrø，2006；侯栋甲等，2014a；杜炳毅等，2016)，使该方法得到了快速发展.由于VTI介质的反演问题引入了三个各向异性参数，导致常规的反演方法往往不能获得精确的各向异性参数反演值，因此VTI介质的反演方法目前只是进行了一些初步的研究且并不完善(刘玉柱等，2014).王明春(2007)利用多波叠前联合反演技术获得了VTI介质的岩性参数和弹性参数；钟峙等(2011)给出了VTI介质的岩性参数反演方法；侯栋甲等(2014b)提出了基于贝叶斯理论的VTI介质多波叠前联合反演方法；刘玉柱等(2015)基于Born敏感核函数实现了一种新的VTI介质多参数全波形反演方法；王义和董良国(2015)采用截断牛顿法在时间域同时反演了VTI介质的声波双参数.传统的多分量叠前联合反演都是将多个目标函数以一定的权重加权叠加构成单目标函数进行优化，从而引入了权重系数对反演结果的影响.此外，该反演问题可以利用多目标函数优化算法有效地解决，这类算法可以在不需要人为地设定权重系数的条件下同时优化多个目标函数，从而得到反演问题的全局最优解，因此具有较强的鲁棒性.目前，许多学者(Kozlovskaya et al., 2007；Singh et al., 2008；Heyburn and Fox, 2010)利用多目标函数优化算法成功地解决了多分量叠前反演问题；Padhi和Mallick(2013, 2014)引入了快速非支配排序遗传算法(NSGA Ⅱ)来处理多分量叠前波形反演问题，成功地获得了地下地层的弹性参数.

本文首先介绍了VTI介质反射率法和NSGA Ⅱ反演方法的基本原理和具体步骤；其次为了提高反演结果的精度，在NSGA Ⅱ算法中引入了线性变换，并且给出了交叉概率、交叉分布指数、变异概率和变异分布指数等遗传算法参数的具体变化形式；最后以反射率法为基础，同时结合NSGA Ⅱ算法实现了VTI介质的多分量叠前联合反演.待求反演参数分别为厚度、纵横波速度、密度和各向异性参数(T，V_P，V_S，ρ，ε和δ)等多个参数.模型测试验证了该反演方法的有效性和正确性.

1 方法原理

对于本文提出的VTI介质多分量叠前联合反演方法来讲，大体可以分为反射率法正演和快速非支配排序遗传算法反演两个部分，因此下面将分别介绍它们的基本理论.为了降低反演的计算量，本文只考虑一维水平层状模型情形并且在截距时间—射线参数域(τ-p)进行整个反演过程.

1.1 VTI介质反射率法正演

在一维水平层状VTI介质模型中，其动量方程和本构方程分别如下：

(1)

(2)

式中，ρ是介质密度；ω是角频率；u是弹性位移矢量；σ是二阶应力张量；e是二阶应变张量；f是单位体积的体力；c是四阶弹性刚度张量.由于对称性，c的81个元素最多有21个是独立的，并且可以表示为6×6的弹性刚度矩阵C.对于VTI介质，其弹性刚度矩阵具体形式为(Auld，1973)

(3)

在右手笛卡尔坐标系下，因为VTI介质的性质仅在垂向上发生变化，所以公式(1)经过傅里叶变换且进行整理后可得到常微分方程：

(4)

式中，A是6×6的系统矩阵，其元素是由弹性常数、水平慢度和密度所构成的函数关系式；F =ω^-1[0, 0, 0, if_x, if_y, if_z]^T；水平慢度

进一步求解常微分方程(4)即可得到VTI介质在频率-慢度域(ω-p)中各个反射界面的位移，然后采用Kennett递推方法(1983)即可计算获得整个地质模型的总反射率.Kennett递推方法的具体形式如下：

(5)

式中，r_d和t_d分别为下行波入射时产生的反射系数、透射系数矩阵；r_u和t_u分别为上行波入射时产生的反射系数、透射系数矩阵；E_u和E_d可称为相位校正项，其具体形式分别为

(6)

式中，λ₁，λ₂和λ₃为系统矩阵A虚部小于零的特征值，对应着上行波；λ₄，λ₅和λ₆为系统矩阵A虚部大于零的特征值，对应着下行波；t_i=ω(z_i-z_i-1)，z_i和z_i-1分别为界面i与i-1的深度.

基于上述所有公式，就可以计算获得ω-p域的总反射系数，然后利用傅里叶反变换将其变换到τ-p域，最后加入震源项就可以获得包含多种地震波传播效应的VTI介质τ-p域全波场地震记录.理论上，反射率法可以同时模拟P波、SV波和SH波，但是在VTI介质中只有P波分量和SV波分量之间存在相互耦合关系，即SH波分量与其他分量无关.因此为了简化研究，本文不考虑SH波分量，即只考虑除γ外的其他五个Thomsen参数.

1.2 快速非支配排序遗传算法反演

多分量叠前联合反演是一类高度非线性问题，理论上它具有非唯一解.目前，解决此类问题一般主要采用目标函数权重加权法进行处理，但是该方法要求人为地设定一组权重系数以构成单目标函数进行优化，因此势必会引入人为因素对反演结果的影响.解决多分量叠前联合反演的另外一种有效方法是多目标函数优化算法，该类方法的主要目的是在满足模型空间或决策空间中一些约束条件的条件下找到一组能够同时最大化或者最小化多个目标函数的折中解，这些解被称为Pareto最优解.快速非支配排序遗传算法就是一种比较高效的多目标函数优化算法(Deb et al., 2002)，该方法是一种带有精英策略的多目标函数全局最小寻优算法，它不需要人为设定权重系数就可以同时最小化多个目标函数，从而得到反演问题的Pareto最优解，最终再根据实际问题的需要从Pareto最优解中提取反演最优解.

为了利用NSGA Ⅱ反演方法从多分量地震记录中反演获得VTI介质的厚度、纵横波速度、密度和各向异性参数等信息，首先应用互相关原理构建如下所示的多目标函数：

(7)

式中，D_PP和D_PSV分别表示实际观测得到的PP波和PSV波地震记录；S_PP和S_PSV分别表示由反演结果正演模拟得到的PP波和PSV波地震记录.

在详细介绍NSGA Ⅱ算法之前，必须首先明白NSGA Ⅱ的非支配排序和种群多样性保护两个主要机制.假定存在一个种群大小为N的种群，首先通过公式(7)可以计算得到种群中每个成员的目标函数值，然后根据该目标函数值对种群中各个成员分配一个非支配等级，该过程即为非支配排序.非支配排序遵循以下准则：非支配等级1的成员支配其它非支配等级的成员，但是它们内部成员之间不互相支配；非支配等级2的成员支配非支配等级3及其以上非支配等级的成员，同样地，同一非支配等级之间也不存在互相支配关系，以此类推直至种群中所有的成员都被分配到一个非支配等级为止.与此同时，为了保持种群的多样性，NSGA Ⅱ为属于同一非支配等级的成员分配一个称为拥挤距离的参数，以此来防止在遗传操作的过程中种群丧失多样性.当种群中每个个体都具有拥挤距离属性时，在后续的锦标赛选择中，同一非支配等级中拥挤距离越大的个体越容易被挑选出来参与到下一次遗传迭代过程中.为了计算拥挤距离，首先需将属于同一非支配等级的成员分别按照多目标函数值进行升序排列，然后对于位于两侧的成员分别分配一个无穷大的数作为其拥挤距离，而位于两侧成员之间的其他成员则可以根据一定的准则计算出不同的拥挤距离.

图 1为NSGA Ⅱ多目标函数反演方法的流程图.从图中可以看出该反演方法的主要步骤可分为：(1)首先在没有违背决策空间或模型空间约束条件下随机生成种群大小为N的初始种群，对其进行快速非支配排序和拥挤距离计算使得种群中每个个体具有非支配等级和拥挤距离两个属性参数；(2)接着利用锦标赛选择、模拟二进制交叉和实参数变异等一系列遗传操作，计算得到种群大小为N的子代种群；(3)然后将父代种群和子代种群合并成为一个种群大小为2N的合并种群，对其进行快速非支配排序和拥挤距离计算，根据其非支配等级和拥挤距离从中挑选出N个新个体作为父代种群参与到下一次遗传操作中；(4)最后重复以上(2)(3)遗传迭代过程直至预先设定的停止准则被满足.

一般地，NSGA Ⅱ算法进行非支配排序和拥挤距离计算主要是根据种群中每个个体的原始目标函数值.然而，由于在遗传迭代的早期种群中绝大多数成员都远离全局最优解，因此纯粹基于原始目标函数值进行这两项工作往往会造成种群中的较优个体驱使NSGA Ⅱ算法向局部最优方向收敛，从而陷入局部极小值(Goldberg，1989).为了避免这一现象的发生，本文采用线性变换对原始目标函数值进行一定的处理，相对降低较优个体的适应度值以及相对升高较劣个体的适应度值，从而使得多目标函数反演过程稳定地向全局最优方向收敛.

线性变换的表达式可以表示为

(8)

和

(9)

式中，v是原始目标函数值；v′是线性变换后的目标函数值；v_avg是原始目标函数值的平均值；v_max是原始目标函数值的最大值；C_m是人为设定的参数，以便保持关系式v′_max=C_mv_avg成立，此时v′_max为线性变换后目标函数值的最大值.为了获得更好的反演结果，在本文的反演过程中设定C_m随着遗传代数逐渐改变，具体形式如下：

(10)

式中，t为当前遗传代数；t_max为最大遗传代数.

传统的遗传算法(Genetic Algorithm，GA)主要采用二进制编码方式对模型参数进行编码，因此只能在模型空间离散采样进行搜索.然而由于反演问题的决策空间或模型空间一般都是连续的，因此为了不丢失全局最优解，本文采用一种在模型空间连续采样的实值编码方式对模型参数进行编码，相对应地引入了Deb和Agrawal(1995, 1999)提出的模拟二进制交叉(Simulated Binary Crossover，SBX)和实参数变异(Real Parameter Mutation，RPM)方法.采用实值编码方式生成的父代种群经过模拟二进制交叉和实参数变异操作后获得子代种群，从而使得NSGA Ⅱ算法能在整个模型空间进行搜索而不丢失全局最优解.

在模拟二进制交叉中，首先定义参数β：

(11)

式中，x_i^{1, t}和x_i^{2, t}分别表示任意一代t中父代个体1和个体2的第i个编码值，一般假定x_i^{2, t}>x_i^{1, t}；x_i^{L, t}和x_i^{U, t}分别表示任意一代t中父代个体模型参数的最小取值和最大取值.

然后定义参数α：

(12)

式中，η为交叉分布指数，它是一个非负实数.若η取值较大，则生成的子代个体靠近其父代个体；若η取值较小，则生成的子代个体远离其父代个体.在理想情况下，遗传迭代开始时η取值应该较小，然后随着遗传代数的增加而逐渐增大.

接着定义参数β_qi：

(13)

式中，u_i是一个介于[0, 1]之间的任意数.

最后可以利用下列公式计算得到两个子代个体：

(14)

式中，x_i^{1, t+1}和x_i^{2, t+1}分别表示任意一代t中子代个体1和个体2的第i个编码值.

在NSGA Ⅱ算法中，为了避免遗传迭代过程陷入局部极小值，经过模拟二进制交叉生成的子代个体必须经过实参数变异操作以跳出局部极小陷阱.在实参数变异中，首先定义参数φ：

(15)

式中，x_i^child表示经过交叉后子代个体的第i个编码值；x_i^L和x_i^U分别表示经过交叉后子代个体模型参数的最小取值和最大取值.

然后定义参数φ_i：

(16)

式中，r_i是一个介于[0, 1]之间的任意数；k为变异分布指数，它是一个非负实数.与η相似，若k取值较大，则生成的变异个体靠近其原始个体；若k取值较小，则生成的变异个体远离其原始个体.在理想情况下，k取值应该随着遗传代数的增加而逐渐改变.

最后可以利用下列公式计算获得经过变异后的子代个体：

(17)

式中，x_i^mutated表示经过变异后子代个体的第i个编码值；x_i^L和x_i^U分别表示经过交叉后子代个体模型参数的最小取值和最大取值.

为了NSGA Ⅱ算法寻优过程的稳定进行，在模拟二进制交叉和实参数变异过程中，应该设定相对应的交叉概率(P_c)和变异概率(P_m)，其取值介于[0, 1]之间，一般地交叉概率取值较大而变异概率取值较小.在本文中，设定这些遗传算法参数(P_c，η，P_m和k)随着代数的增加而逐渐改变(Li and Mallick, 2015)：

(18)

式中，t为当前遗传代数；t_max为最大遗传代数；n为总变量数.

2 模型测试

为了说明本文VTI介质多分量叠前联合反演方法的正确性、有效性以及可行性，本文采用两个不同的一维VTI介质模型测试该反演方法的性能.每次测试过程中都考虑了不同因素对反演结果的影响，并对其反演结果进行了对比分析.

2.1 简单理论模型

为了验证本文反演方法的正确性和有效性，本文首先设计了一个简单的层状VTI介质理论模型(模型A)进行测试.该模型分为7层，总共440 m，每层的厚度均不相同，最大层厚和最小层厚分别为120 m和5 m.表 1为该模型的各层参数.本文采用反射率法计算该模型的PP波和PSV波合成地震记录，子波采用主频为30 Hz的雷克子波，射线参数p取值介于0~0.2 s·km^-1之间，间隔为0.01 s·km^-1. 图 2为正演模拟得到的未含噪声多分量合成地震记录，图 2a为PP波地震记录，图 2b为PSV波地震记录.考虑到反射率法在正演过程中考虑了地层厚度信息，所以本文针对该模型将地层厚度也作为待求反演参数，因此其反演结果分别为地层厚度、纵横波速度、密度和各向异性参数6个参数，则需要反演的模型参数矢量共由42个未知数组成.在该模型的反演过程中，本文设定种群大小为800，最大遗传迭代次数为1000.

针对如图 2所示的多分量合成地震记录，本文研究了以下三方面内容：(1)利用传统的单目标函数反演方法和本文的多目标函数反演方法分别进行反演，对比分析了它们的反演结果；(2)为了分析噪声对该多分量叠前联合反演方法的影响，本文对如图 2所示的多分量地震记录分别加入了不同信噪比(SNR=2，4，6，8)的随机噪声进行反演，对比分析了它们的反演结果；(3)利用不同的搜索窗对反演过程进行约束，讨论了搜索窗对反演结果的影响.图 3为使用NSGA Ⅱ反演方法时模型A的厚度搜索窗(搜索范围)，图 3a对应于优化搜索窗，图 3b对应于垂直搜索窗.

图 4为基于传统单目标函数反演方法的各参数反演结果，黑色实线为初始模型，蓝色实线为真实值，红色实线为反演值.在该反演过程中，本文采用Graebner(1992)推导的VTI介质反射系数精确公式作为正演方程，由于该公式在正演的过程中未考虑地层厚度，因此首先需要将模型A转换到时间域进行离散采样，设定时间采样间隔为2 ms；此外还需将如图 2所示的多分量地震记录进行纵横波匹配和同相轴拉平处理.由图 4可知，基于传统单目标函数反演方法的各参数反演结果在整体上接近于真实值，但是由于该方法未考虑地层厚度，所以不同时间采样点的反演结果会在真实值附近出现不同偏差，表现为同一层的反演结果在真实值附近出现“跳动”.当加入信噪比为2的随机噪声时，其反演结果的误差进一步增大，即反演结果的精度降低.另外从图中还可以看出，该传统单目标函数反演方法对薄层的识别能力不足.图 5和图 6分别为不同信噪比的情形下使用优化搜索窗和垂直搜索窗时基于NSGA Ⅱ反演方法的模型A各参数的反演结果.

从图 5可以看出，对于层状地质模型来讲，由于本文反演方法在反演过程中考虑了地层厚度，因此其反演结果未出现如传统单目标函数反演方法反演结果的“锯齿状”现象.对比分析图 5a—5e可知，不同信噪比的反演结果在整体上都趋近于模型A的真实值.多分量地震数据未含噪声时，模型参数的反演结果与真实值之间的误差最小，薄层也能被准确地识别，反演结果的精度最高，但随着数据信噪比的逐渐降低，模型参数的反演值与真实值之间的误差逐渐增大，对于薄层的识别能力也逐渐变弱，即反演结果的精度逐渐降低.从图 6可以看出，与图 5类似，当多分量地震数据未含噪声时，各参数的反演结果最为准确，然而随着数据信噪比的逐渐降低，模型参数反演结果的精度也会逐渐降低.从图 6e可以看出，当地震数据的信噪比为2且使用垂直搜索窗时，各参数的反演结果仅仅在变化趋势上与模型参数真实值基本保持一致，反演结果与真实值之间误差较大，薄层反演结果的精度也较低.对比分析图 5和图 6可以得出，基于优化搜索窗的反演结果在整体上优于基于垂直搜索窗的反演结果，这一现象对于薄层更为明显，即相对于垂直搜索窗，使用优化搜索窗能更为准确地反演获得薄层信息，这是因为优化搜索窗对于多分量叠前联合反演过程的约束能力强于垂直搜索窗，尤其是对薄层.

为了进一步分析该多分量叠前联合反演方法的性能，本文对模型A各参数的反演结果进行了误差分析，建立了如表 2—6所示的五个反演结果的误差统计表.表 2为基于传统单目标函数反演方法的各参数反演结果的平均误差统计表.从表中可以看出，当地震数据未含噪声时，模型A纵横波速度和密度反演结果的相对误差均未超过8.26%，但两个各向异性参数的相对误差分别为21.90%和24.41%，说明该方法的纵横波速度和密度反演结果的精度较高，各向异性参数反演结果的精度较低.当加入信噪比为2的随机噪声时，各参数反演结果的相对误差增大，此时其纵横波速度和密度反演结果的相对误差均未超过12.33%，两个各向异性参数的相对误差分别增大到47.68%和49.81%，说明噪声对各向异性参数影响较大.表 3为不同信噪比的情形下使用优化搜索窗时模型A各层的厚度反演结果和绝对误差统计表.从表中可以看出，当多分量地震数据未含噪声时，模型A的各层厚度反演结果与真实值之间的绝对误差整体基本最小，即与真实值最为接近.随着数据信噪比的逐渐降低，模型A的各层厚度反演结果与真实值之间的绝对误差整体上逐渐增大，即厚度反演结果的精度整体上逐渐降低.表 4为不同信噪比的情形下使用优化搜索窗时模型A的各参数反演结果的绝对误差和相对误差统计表.从表中可以看出，随着数据信噪比的逐渐提高，模型A各参数的反演结果与真实值之间的绝对误差和相对误差整体上都呈现逐渐减小的趋势，当数据未含噪声时，反演结果的绝对误差和相对误差基本达到最小.当使用优化搜索窗约束不同信噪比的反演过程时，模型A纵横波速度和密度反演结果的相对误差均未超过5.31%，说明这三个参数的反演精度较高；模型A厚度和各向异性参数反演结果的相对误差随着信噪比的逐渐降低，其值急剧增大，但都未超过27.67%，结合其绝对误差分析，它们的反演结果是合理的.综上，说明噪声对纵横波速度和密度的影响较小，对厚度和各向异性参数的影响较大.表 5为不同信噪比的情形下使用垂直搜索窗时模型A各层的厚度反演结果和绝对误差统计表.从表中可以看出，与表 3类似，对于未含噪声的多分量地震数据，其反演结果与真实值之间的绝对误差也整体基本最小，表明其厚度反演结果符合真实值.随着数据信噪比的逐渐降低，厚度反演结果与真实值之间的绝对误差也整体上逐渐增大.对比分析表 3和表 5可得出，基于优化搜索窗的厚度反演结果整体上优于基于垂直搜索窗的厚度反演结果.表 6为不同信噪比的情形下使用垂直搜索窗时模型A的各参数反演结果的绝对误差和相对误差统计表.从表中可以看出，与表 4类似，当多分量地震数据未含噪声时，各参数反演结果的绝对误差和相对误差也基本达到最小.随着数据信噪比的逐渐降低，各参数的反演结果与真实值之间的绝对误差和相对误差也逐渐增大.当使用垂直搜索窗约束不同信噪比的反演过程时，模型A纵横波速度和密度反演结果的相对误差均未超过8.80%，同样说明这三个参数的反演精度较高；模型A厚度和各向异性参数反演结果的相对误差随着信噪比的逐渐降低，其值同样急剧增大，但都未超过42.92%，结合其绝对误差分析，它们的反演结果也是合理的.综上，同样说明噪声对纵横波速度和密度的影响较小，对厚度和各向异性参数的影响较大.对比分析表 2、表 4和表 6可以得出，相对于传统的单目标函数反演方法，通过基于NSGA Ⅱ算法的多分量叠前联合反演方法反演得到的各参数反演结果的相对误差在整体上要小一些，尤其是各向异性参数，这说明本文的反演方法能更为准确地反演获得VTI介质的各向异性参数.基于优化搜索窗的各参数反演结果在整体上优于基于垂直搜索窗的各参数反演结果，说明一个好的搜索窗可以提高反演结果的精度.

图 7为不同情形下模型A在反演过程中的目标函数变化曲线.图 7a和图 7c分别对应于当多分量地震数据未含噪声时使用优化搜索窗和垂直搜索窗约束反演过程的PP波和PSV波收敛曲线，图 7b和图 7d分别对应于当多分量地震数据含有噪声(SNR=2)时使用优化搜索窗和垂直搜索窗约束反演过程的PP波和PSV波收敛曲线.从图中可以看出，在前100次遗传迭代过程中，该联合反演方法的收敛速度较快，随后，随着遗传代数的增加，其收敛速度逐渐降低并趋于平缓.对于无噪声地震记录，其PP波和PSV波的收敛程度最好，当数据含有噪声时，PP波和PSV波的收敛程度变差，这一现象在PSV波上表现更为明显.此外对比分析图 7a—7d可以发现，使用优化搜索窗时PP波和PSV波的收敛程度略好于垂直搜索窗.图 8为不同情形下模型A在反演过程中的最后一代在目标空间的种群分布图，实际上它们反映了通过NSGA Ⅱ算法全局寻优最终得到的Pareto最优前沿.图 8a和图 8c分别对应于当多分量地震数据未含噪声时使用优化搜索窗和垂直搜索窗约束反演过程得到的Pareto最优前沿，图 8b和图 8d分别对应于当多分量地震数据含有噪声(SNR=2)时使用优化搜索窗和垂直搜索窗约束反演过程得到的Pareto最优前沿.从图中可以看出，当地震数据未含噪声时，Pareto最优前沿的收敛程度最好，即最接近于全局最优解，当地震记录含有噪声时，其收敛程度会随之降低，同时Pareto最优前沿收敛性的变差程度在PSV波上更为剧烈.此外同样可以发现使用优化搜索窗时Pareto最优前沿的收敛程度略优于垂直搜索窗.综上可以得出，噪声会降低该多分量叠前联合反演方法的收敛程度，尤其是对PSV波，当多分量地震数据没有噪声时，PP波和PSV波的收敛程度最好.同时，一个好的搜索窗可以在一定程度上提高该反演方法的收敛程度，从而得到更优的反演结果.

2.2 一维多层模型

为了进一步验证本文反演方法的可行性，本文基于某实际测井数据建立了一个一维多层VTI介质模型(模型B)对该叠前联合反演方法进行测试.该模型共有50层，每层的厚度设定为50 m，则该模型总共2500 m.图 9显示了该VTI介质模型的各层参数.图 10为利用反射率法计算的该模型在τ-p域的PP波和PSV波合成地震记录，其中雷克子波主频为30 Hz，射线参数p取值介于0~0.15 s·km^-1之间，间隔为0.01 s·km^-1.在此模型的反演过程中，本文设定种群大小为800，最大遗传迭代次数为2000.

针对如图 10所示的PP波和PSV波合成地震记录，本文主要研究了以下三方面内容：(1)分别设定地层厚度已知和未知两种情况进行反演，对应的待搜索模型参数矢量分别含有250和300个未知数，对比分析了它们的反演结果；(2)通过加入信噪比为2(SNR=2)的随机噪声进行反演，分析讨论了噪声对反演结果的影响；(3)最后同样分析讨论了搜索窗对反演结果的影响.对于模型B，本文使用了两个搜索窗来约束其反演过程，其中之一(搜索窗A)由模型参数的真实值向两侧展宽得到，此时设定厚度搜索范围为40~60 m；另一个(搜索窗B)则是通过模型参数的真实值平滑10次后向两侧展宽得到，此时设定厚度搜索范围为30~70 m，显然搜索窗A优于搜索窗B.

图 11是未含噪声时使用搜索窗A约束反演过程得到的模型B的反演结果，其中图 11a和图 11b分别显示了其已知地层厚度和未知地层厚度两种情形下的反演结果曲线与真实值曲线.从图中可以看出，当地震数据未含噪声时，反演结果曲线与真实值曲线整体吻合较好.当同时反演地层厚度时，其厚度反演结果误差最大不超过6.1 m，由于厚度反演结果存在误差，使得某些深度处其他参数的反演结果超出了其对应的搜索范围；另外，与已知地层厚度相比，虽然各参数反演结果的整体变化趋势与真实值仍然大体一致，但各参数的反演值与真实值之间的误差增大.图 12是含噪声(SNR=2)时使用搜索窗A约束反演过程得到的模型B的反演结果，其中图 12a和图 12b分别显示了其已知地层厚度和未知地层厚度两种情形下的反演结果曲线与真实值曲线.同样可以从图中看出，当地层厚度未知时各参数反演结果的精度会有所下降，此时厚度反演结果误差最大不超过8.4 m.对比分析图 11和图 12可知，当信噪比为2时，各参数反演结果的精度降低，尤其各向异性参数更为明显，说明噪声对各向异性参数的影响较大.图 13和图 14分别类似于图 11和图 12，不同的是它们分别对应于使用搜索窗B时未含噪声和含噪声(SNR=2)的模型B反演结果，其中当厚度未知时其厚度反演结果的最大误差分别不超过10.6 m和16.7 m.类似地，可以从图中看出，地层厚度未知和噪声会在一定程度上降低反演结果的精度，噪声对厚度和各向异性参数影响较大.对比分析图 11—14可以发现，搜索窗A的反演结果整体上优于搜索窗B的反演结果，厚度搜索窗对厚度反演结果影响较大.值得注意的是，本文设定模型B各层厚度均为50 m，因此该反演方法能够获得较为准确的厚度反演值；若地层厚度逐渐变薄，其厚度反演结果的精度势必会随之降低.对于全为薄层的地质模型，该反演方法恢复其各薄层厚度信息的能力仍然不足，因此对于该类情况一般假定地层厚度已知进行反演，通常设定为某一恒定的深度采样间隔.

图 15为在地层厚度未知且未含噪声的情形下使用搜索窗A约束模型B反演过程时种群随遗传代数增加的演化，其中横纵坐标分别代表着如公式(7)所定义的PP波和PSV波目标函数值经过线性变换后的值.从图中可以看出，随着遗传代数的增加，种群中各个成员的PP波和PSV波误差同时逐渐靠近坐标原点，即种群逐渐向误差最小方向收敛，说明该多目标函数反演方法是有效的.另外可以发现，该反演过程在前500次的迭代过程中其收敛速度较快，随后随着代数的增加其收敛速度逐渐变慢直至趋于平缓.图 16为未含噪声的情形下使用搜索窗A时模型B在反演过程中的最后一代在目标空间的种群分布，图 16a和16b分别对应于已知地层厚度和未知地层厚度两种情形.从图中可以看出，当地层厚度未知时，即同时反演地层厚度时，由该多目标函数反演方法得到的Pareto最优前沿的收敛程度差于已知地层厚度的情形，说明地层厚度未知时各参数反演结果的精度会降低.对比分析模型A和模型B的反演过程可以发现，当模型的层数增多时，即当待反演模型参数矢量的未知数增多时，NSGA Ⅱ多目标函数反演方法的收敛速度和收敛程度同样都会变差.

3 结论

本文通过反射率法和快速非支配排序遗传算法研究了VTI介质在τ-p域的多分量叠前联合反演方法.为了能够得到更好的反演结果，本文将线性变换引入到NSGA Ⅱ反演方法中对种群中各个成员的原始目标函数值进行适当的调整，同时设定在反演的过程中逐渐改变各个遗传算法参数，促使整个反演过程稳定地向全局最优方向收敛.

理论模型的反演测试验证了本文反演方法的正确性、有效性和可行性.与传统的多分量叠前联合反演方法相比，本文的反演方法不仅能够获得较为精确的VTI介质纵横波速度、密度和各向异性参数，同时还可以反演得到地层的厚度信息；但是如果将地层厚度同时作为待求的反演参数，则其他参数的反演精度就会有所降低，同时该反演方法对薄层的识别能力受噪声和搜索窗的影响较大.不同信噪比的多分量地震记录的反演结果说明了本文的反演方法具有良好的抗噪性，噪声对厚度和各向异性参数的影响较大.不同搜索窗的反演结果表明，搜索窗在一定程度上影响着反演结果的精度，随着搜索窗取值逐渐偏离真实模型，反演结果的精度会有所下降.随着反演未知数或目标函数的增加，本文基于NSGA Ⅱ算法的多分量叠前联合反演方法的收敛性会变差以及其计算量也会急剧增加.

VTI介质的多分量叠前联合反演是一项极其复杂的研究工作，需要注意的方面有很多.例如，本文只是实现了一维模型的反演测试，有关更为复杂的模型反演有待进一步研究；本文只进行了两种简单搜索窗的对比分析，但在实际的反演过程中，如何利用已获的地质、地震以及测井等资料建立搜索窗同样需要进一步研究；相对于其他非种群进化类的单目标函数优化方法，该反演方法的计算量仍然较大，因此如何降低其计算量也是一项重要的研究工作.

附录A

包含NSGA Ⅱ在内的所有多目标函数演化算法(Multi-objective Evolutionary Algorithms，MOEA)的主要目的都是同时最大化或者最小化多个目标函数从而获得多目标优化问题的Pareto最优解集.为了更为完整地阐述本文VTI介质的多分量叠前联合反演方法，下面将详细介绍NSGA Ⅱ算法的一些主要概念.

1 Pareto最优解、Pareto最优解集以及Pareto最优前沿

多分量叠前联合反演一般可以归纳为解决多目标函数优化问题(Multi-objective Optimization Problem，MOP).通常多目标优化问题的最终目的是在一组约束条件下，同时最大化或者最小化多个不同的目标函数，其数学表达形式为

(A1)

式中，X=[x₁, x₂, …, x_l]是一个l维向量；f_i(X), i=1, 2, …, n是目标函数；g_j(X)和h_k(X)为约束条件.

一般地，多目标优化问题与单目标优化问题之间存在着很大差异.对于单目标优化问题，因为其只有一个目标函数，因此通常需要寻找优于其他所有解的最好解，一般为全局最大或者全局最小.而当存在多个目标函数时，由于目标函数之间存在冲突或者无法比较的现象，所以往往很难找到一个解使得所有的目标函数同时达到最小或者最大，也就是说，一个解可能对于某个目标函数是最好的，但是对于其他目标函数可能却不是最好的，甚至是最差的.因此对于多目标优化问题，通常存在一系列折中解，就当前的目标空间而言这些解是无法比较它们优劣的，因此这种解被称为非支配解(Nondominated solutions)或Pareto最优解(Pareto-optimal solutions).如果用x^*进行表示，则需满足以下条件：

(a) F(x)=min[f₁(x), f₂(x), …, f_n(x)]；

(b) F(x^*)=min[f₁(x^*), f₂(x^*, …, f_n(x^*];

(c) 不存在x ∈ Ω，使得F(x) < F(x^*).

式中，Ω为变量空间或决策空间；F为对应的目标空间.

所有Pareto最优解组成的集合被称为Pareto最优解集(Pareto-optimal solution set).Pareto最优前沿(Pareto-optimal front)则为Pareto最优解集在对应的目标空间中的像.

2 支配与非支配

对于拥有m个目标函数的多目标函数优化问题，如果一组解x¹支配另一组解x²，就可以表示为

(A2)

因为NSGA Ⅱ是一种最小化寻优算法，因此当满足以下条件时，就可以认为如果当决策空间的两组解x¹=[x₁¹, x₂¹, …, x_n¹]^T和x²=[x₁², x₂², …, x_n²]^T分别映射到目标空间对应着y¹=[y₁¹, y₂¹, …, y_m¹]^T和y²=[y₁², y₂², …, y_m²]^T时，且满足

(a) y_i¹≤y_i²，i=1, 2, …, m；

(b)至少存在一个j使得y_j¹ < y_j²，其中1≤j≤m.

3 拥挤距离

拥挤距离的主要目的是为了避免某一较优个体在寻优迭代的过程中快速充满整个种群导致寻优过程过早收敛而陷入局部极小值.最初，非支配排序遗传算法(NSGA)通过人为地设定共享参数σ_share来达到这一目的，而快速非支配排序遗传算法(NSGA Ⅱ)首次提出了拥挤距离.拥挤距离是指在种群中某成员的周围个体的密度，可以用i_d表示，也可以简单理解为在个体i周围只包含个体i本身的最大长方形的长，如图A1.

4 锦标赛选择

当种群完成非支配排序和拥挤距离计算时，种群中的每个个体都具备了非支配等级i_rank和拥挤距离i_d两个属性参数.为了使寻优问题向最优方向进化发展，需要利用锦标赛选择方法根据这两个属性从该种群中选择父代种群进行后续的交叉变异操作.在锦标赛选择的过程中，首先需要从种群中随机挑选多个成员，然后对挑选出来的成员进行优劣对比，只保留其中最优的个体存入父代种群，其他的个体则被舍去.其基本原则是当个体i优于个体j时需至少满足下列条件的其中之一：

(a) i_rank < j_rank；

(b) i_rank=j_rank且i_d>j_d.

5 精英策略

为了防止在遗传迭代的过程中丢失最优解，NSGA Ⅱ算法提出了精英保护机制，即精英策略.当子代种群生成时，首先将种群大小均为N的父代种群和子代种群合并成为一个种群大小为2N的合并种群；然后对该合并种群进行非支配排序和拥挤距离计算；最后从该合并种群中挑选出N个新个体作为父代种群参与到下一次遗传操作中.其挑选的基本原则为：首先属于非支配等级1的所有成员均被选择，然后依次类推直至加入某一非支配等级的所有成员时新父代种群的大小大于N时停止；此时为了充满新父代种群，首先将该非支配等级的所有成员按照拥挤距离进行降序排列，然后根据拥挤距离的大小从大到小依次挑选出较优成员加入到新父代种群直至其种群大小等于N.

通过该精英策略操作，可以保留那些较优个体，舍去那些较劣个体，从而使得多目标优化问题向全局最优方向收敛.

致谢  感谢中国石油大学(北京)李景叶教授对本项研究工作的支持与帮助！感谢两位审稿专家对本文提出的指导和建议！