CN102759746A - 一种变偏移距垂直地震剖面数据反演各向异性参数方法 - Google Patents

Abstract

本发明是地震勘探的变偏移距垂直地震剖面数据反演各向异性参数方法。计算零偏垂直地震剖面各层的垂向纵、横波速度，拾取标志层位确定各层与井轨迹交点，建立表层模型和观测系统，随机初始化待求层的各向异性参数和待求层顶界面的3次系数，形成待优化的各向异性速度模型，筛选层之间的检波点，拾取炮点初至计算极化角，使用P波初至旅行时和极化角建立目标函数，利用小生境遗传算法优化目标函数，完成待求层各向异性参数和待求层顶界面形态的优化，“从浅到深剥层”逐层计算完成整个模型各向异性参数优化。本发明反射界面更接近真实地层，采用线性收敛的射线追踪计算效率高，简化了各向异性介质复杂的内在关系，小生境遗传算法导向搜索指导下迭代试射射线追踪正演，非线性反演优化模型各向异性参数、界面形态。

Description

一种变偏移距垂直地震剖面数据反演各向异性参数方法

技术领域

本发明涉及地震勘探数据处理技术，是一种变偏移距垂直地震剖面(Walkaway VSP)数据反演各向异性参数方法。

背景技术

Walkaway VSP采用井中接收，震点在地面沿一条过井线移动的观测方式，主要用于井旁构造成像及储层研究。3D VSP(3维垂直地震剖面)可以视为不同方位的Walkaway VSP组成。Walkaway VSP数据包含着丰富的速度信息和反射界面形态参数，表现为波至旅行时和极化角的变化。

Thomsen(1986)提出了具有明确物理意义的各向异性参数，给出了任意强度和弱各向异性的相速度和群速度表达式，使各向异性介质研究进入了应用阶段。Gelinsky等(1999)用二维射线追踪计算VTI(具有垂直对称轴的横向各向同性)介质模型Walkaway VSP初至旅行时，精细调整模型速度参数，反复迭代直至模型初至旅行时与实际初至旅行时拟合差最小，得到二维平层VTI模型。Song等(2001)提出利用联合旅行时和极化数据的弹性波反演方法，极化数据应用的线性化问题从射线扰动理论着手，并假定介质是弱各向异性。Zhou等(2003)在WalkawayVSP处理中联合应用局部慢度和旅行时反演，利用射线追踪比较测量的和各向异性速度模型计算的旅行时，通过先导试验和误差技术给出测量旅行时的最小平方拟合。Zhou等(2004)用有限差分算法计算二维TTI介质Walkaway VSP初至旅行时，迭代正演估算TTI介质的等效Thomsen参数，计算精度用旅行时均方差表示。Zhou等(2006)采用Walkaway VSP旅行时线性反演确定Thomsen各向异性参数ε和δ，反演二维垂直对称轴的横向各向同性(VTI)，模型旅行时与观测旅行时差为最小平方意义下最小值。

罗小明等(2005)联合井地资料求取地层各向异性参数。通过零偏VSP井的下行初至波拾取和速度分析，得到火成岩层垂直速度V_P0；通过地表地震近偏移距资料各向同性速度分析获得火成岩的动校正速度V_NMO，利用动校正速度与垂直速度的关系求取δ；通过偏移距VSP井下行初至波速度分析求得相速度，由空间几何关系计算相角，利用Thomsen给出的相速度公式计算ε；进而求得瞬时各向异性参数η。

综上，利用初至旅行时和极化数据反演各向异性模型是可行的。但是，用二维平层VTI介质或等效二维起伏TTI介质不能准确地表示地下介质；线性近似不能准确地表示各向异性模型参数与初至旅行时和极化数据的关系，建立在线性近似基础上的线性反演也不能有效恢复各向异性模型参数与初至旅行时和极化数据的非线性关系；基于弱各向异性假设的方法不具有普适性；线性反演法误差大，枚举法搜索能力有限。

发明内容

本发明目的是提供一种克服弱各向异性介质假设的不足，二维平层VTI介质或二维等效TTI介质表示地层和线性反演(线性近似)的局限性，利用变偏移距垂直地震剖面数据反演各向异性参数的方法。

本发明采用以下技术步骤：

1)由变偏移距垂直地震剖面野外观测系统确定以井口为中心的模型范围，假定条件是各向异性模型为横向各向同性(VTI)介质；

2)用零偏垂直地震剖面(VSP)纵、横波初至旅行时计算垂向纵、横波速度，根据纵波速度的变化趋势划分地层并给出各层与井轨迹的交点，计算各层的层间垂向纵、横波速度；

步骤2)计算垂向横波速度时，若无零偏横波资料，则利用转换横波资料。

3)在变偏移距垂直地震剖面炮线方向的深度域过井地震剖面上，拾取标志层的层位，计算各层与井轨迹的交点；由距离最近的上下两个标志层，采用线性加权计算得到中间层界面的坐标(x_i，z_i)，由这些坐标拟合出界面的3次多项式函数；

步骤3)所述的线性加权函数是：

$z_i=\frac{(z^i-z_m^j)}{(z_m^{j+1}-z_m^j)}{z}_{m,j}+\frac{(z_m^{j+1}-z^i)}{(z_m^{j+1}-z_m^j)}{z}_{m,j+1}---\left(1\right)$

x_i＝x_m，j

式中：(x_i，z_i)是中间层i的坐标；i是地层号；j是标志层号；(xⁱ，zⁱ)是层i与井轨迹交点的坐标；

是标志层与井轨迹交点的纵坐标；

z_m，j、z_m，j+1是拾取的标志层纵坐标；

x_m，j是拾取的标志层的横坐标；

由中间层i的坐标(x_i，z_i)拟合出界面的3次多项式函数为：

z＝a₁·x³+a₂·x²+a₃·x+a₄    (2)

式中：a₁、a₂、a₃、a₄是界面的3次系数；

4)由已知的表层速度资料建立表层模型，根据变偏移距垂直地震剖面野外观测情况给模型建立相应的观测系统；

步骤4)若无表层速度资料，则利用已知的静校正数据，建立不带表层结构的从某一基准面开始的模型，并根据变偏移距垂直地震剖面野外观测情况给模型建立相应的观测系统；

5)在参数的范围内，随机初始化待求层的各向异性参数和待求层顶界面的3次系数，形成待优化的各向异性速度模型；

6)根据待求层顶界面、底界面与井轨迹的交点，筛选位于该层之间的检波点，拾取这些检波点接收到得各个炮点的初至、计算极化角；

7)试射射线追踪正演步骤5)模型在步骤6)的检波点位置接收到的各个炮点初至、极化角；

步骤7)所述的极化角是射线到达检波点时，射线方向与垂向所成的角度。

步骤7)所述的试射射线追踪正演初至旅行时、极化角采用以下方法：

(1)给定初始相角；

(2)计算相速度、群速度、群角；

由相角计算相速度的公式为：

$V_P\left(\mathrm{\θ}\right)=V_{P0}\sqrt{1+\mathrm{\ϵ}{\sin }^2\mathrm{\θ}+\frac{1}{2}(1-\frac{V_{S0}^2}{V_{P0}^2})\left\{\right[1+\frac{4{\mathrm{\δ}}^{*}}{(1-V_{S0}^2/V_{P0}^2)}s{\mathrm{in}}^2\mathrm{\θ}{\cos }^2\mathrm{\θ}+\frac{4(1-V_{S0}^2/V_{P0}^2+\mathrm{\ϵ})\mathrm{\ϵ}}{(1-V_{S0}^2/V_{P0}^2)}{\sin }^4\mathrm{\θ}{]}^{1/2}-1\}}---\left(3\right)$

式中：θ是相角，ε、δ^*是Thomsen各向异性参数，V_P0、V_S0是纵、横波垂直各向同性面的相速度，V_P(θ)是相角θ对应的纵波相速度。

由相角计算群速度的公式为：

$V_g\left(\mathrm{\φ}\right)=\sqrt{V^2\left(\mathrm{\θ}\right)+{\left(\frac{\mathrm{dV}\left(\mathrm{\θ}\right)}{\mathrm{d\θ}}\right)}^2}---\left(4\right)$

式中：θ是相角，φ是群角，V(θ)是相速度，V_g(φ)是群速度。

由相速度、群速度计算群角的公式为：

$\mathrm{\φ}=\arccos \left(\frac{V_g\left(\mathrm{\φ}\right)}{V\left(\mathrm{\θ}\right)}\right)+\mathrm{\θ}---\left(5\right)$

式中：θ是相角，φ是群角，V(θ)是相速度，V_g(φ)是群速度。

(3)射线沿群角方向传播，计算射线与界面的交点；

(4)在界面交点处由斯奈尔定律做透射处理，计算透射相角；

所述的斯奈尔定律为：

$P=\frac{\sin \left({\mathrm{\α}}_{\mathrm{inc}}\right)}{V\left({\mathrm{\θ}}_{\mathrm{inc}}\right)}=\frac{\sin \left({\mathrm{\α}}_{\mathrm{tra}}\right)}{V\left({\mathrm{\θ}}_{\mathrm{tra}}\right)}---\left(6\right)$

式中：α_dip是界面倾角；θ_inc是入射相角，α_inc是入射角；θ_tra是透射相角，α_tra是透射角；

界面交点处的各个角度关系为：

α_inc＝|α_dip+θ_inc|

θ_tra＝sign(α_dip+θ_inc)α_tra-α_dip    (7)

式中：α_dip是界面倾角；θ_inc是入射相角，α_inc是入射角；θ_tra是透射相角，α_tra是透射角；

(5)重复(2)-(4)逐层计算直到射线到达井轨迹；

(6)计算射线与井轨迹的交点，计算交点与检波点的距离，若交点与检波点的距离小于给定精度，该射线即为给定精度条件下的准确射线，计算极化角；

若交点与检波点的距离大于给定精度，则由检波点附近的两条射线线性插值计算初始相角，重复(2)-(6)；

线性插值计算初始相角的公式为：

$\mathrm{\α}={\mathrm{\α}}_1+\frac{z-z_1}{z_2-z_1}({\mathrm{\α}}_2-{\mathrm{\α}}_1)---\left(8\right)$

式中：α₁、α₂是检波点附近的两条射线的初始相角；

(x₁，z₁)、(x₂，z₂)是α₁、α₂对应出射坐标(即射线与井轨迹的交点)；

(x，z)是检波点坐标，α是待求的初始相角。

(7)计算射线路径对应的旅行时；

由射线路径计算旅行时的公式为：

$T_{\mathrm{tra}}=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{l_i}{V_{\mathrm{gi}}}---\left(9\right)$

式中，n是射线段数；

l_i是射线追踪到的第i个射线段；

V_gi是第i个射线段对应的群速度；

T_tra实现路径对应的旅行时；

8)采用P波初至旅行时或P波极化角或联合使用P波初至旅行时和P波极化角建立目标函数；

步骤8)所述的采用P波初至旅行时建立的目标函数为：

$\mathrm{objv}=\sqrt{\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(T_{\mathrm{for}}^k-T_{\mathrm{med}}^k)}^2}{n}}---\left(10\right)$

式中：

n是实测的初至个数；

是该层检波点野外实测的第k个初至；

是该层检波点接收到的模型正演的第k个初至。

步骤8)所述的采用P波极化角建立的目标函数为：

$\mathrm{objv}=\sqrt{\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{\α}}_{\mathrm{for}}^k-{\mathrm{\α}}_{\mathrm{med}}^k)}^2}{n}}---\left(11\right)$

式中：

n是实测的初至个数；

是该层检波点野外实测的第k个极化角；

是该层检波点接收到模型正演的第k个极化角。

步骤8)所述的联合使用P波初至旅行时和P波极化角建立的目标函数为：

$\mathrm{objv}=C_T\sqrt{\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{(T_{\mathrm{for}}^k-T_{\mathrm{med}}^k)}^2}{n}}+C_{\mathrm{\α}}\sqrt{\frac{\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{({\mathrm{\α}}_{\mathrm{for}}^k-{\mathrm{\α}}_{\mathrm{med}}^k)}^2}{n}}---\left(12\right)$

式中：

n是实测的初至个数；

该层检波点野外实测的第k个初至；

该层检波点接收到模型正演的第k个初至；

该层检波点野外实测的第k个极化角；

该层检波点接收到模型正演的第k个极化角；

C_T、C_α是初至、极化角在目标函数中的权系数。

9)利用小生境遗传算法优化步骤8)的目标函数；

步骤8)目标函数最小值对应的步骤5)的待优化的各向异性速度模型就是最优模型，这就得到了待求层的层各向异性参数和该层顶界面的3次系数。

步骤9)中所述的小生境遗传算法的参数是：

编码：实数编码；

初始种群生成：在变量范围，随机生成；

选择算子：随机遍历采样选择；

交叉操作：均匀算术交叉，自适应交叉概率；

变异操作：实值变异，自适应变异概率；

群体间的迭代策略-适值共享排挤：将父代+子代形成新的种群，适值共享调整个体的适值，然后执行排挤操作直至种群个体数与父代相同；

终止条件：最大遗传代数+目标函数小于给定精度。

步骤9)中所述的小生境遗传算法的实现步骤为：

(1)随机产生初始化群体，对个体适应度进行评价，

个体适应度评价的函数为：

$\mathrm{FitnV}=2-\mathrm{sp}+2\×(\mathrm{sp}-1)\×\frac{\mathrm{Pos}-1}{\mathrm{Nind}-1}---\left(13\right)$

式中，sp是压差数；

Pos个体在种群中降序排序的位置；

Nind种群中个体的个数；

FitnV个体适应度；

(2)随机遍历采样选择配对父代个体，

随机遍历采样选择的概率为：

$F\left(x_i\right)=\frac{f\left(x_i\right)}{\underset{i=1}{\overset{\mathrm{Nind}}{\mathrm{\Σ}}}f\left(x_i\right)}---\left(14\right)$

式中，f(x_i)是个体x_i的适应度；

Nind种群中个体的个数；

F(x_i)是个体x_i被选择的概率；

(3)按照自适应概率，均匀算术交叉选择配对的父代个体，

自适应交叉概率为：

$P_c=\left\{\begin{array}{cc}\frac{k_1(f_{\max }-f^{\&prime;})}{f_{\max }-\stackrel{\&OverBar;}{f}}& f^{\&prime;}\&GreaterEqual;\stackrel{\&OverBar;}{f}\\ k_3& f^{\&prime;}<\stackrel{\&OverBar;}{f}\end{array}\right.---\left(15\right)$

式中：f_max是群体中最大的适应值，

是群体的平均适应值，f′是交叉的两个个体中较大的适应值，k₁，k₃是(0，1)间的常数，P_c是交叉概率；

均匀算术交叉函数为：

$X_A^{t+1}=\mathrm{\&alpha;}{X}_B^t+(1-\mathrm{\&alpha;})X_A^t$ (16)

$X_B^{t+1}=\mathrm{\&alpha;}{X}_A^t+(1-\mathrm{\&alpha;})X_B^t$

式中：

是待交叉的父代个体；

α是[0，1]之间的一个常数；

是交叉后的子代个体；

(4)按照自适应概率，实值变异交叉过的个体，产生新的种群，

自适应变异概率为：

$P_m=\left\{\begin{array}{cc}\frac{k_2(f_{\max }-f)}{f_{\max }-\stackrel{\&OverBar;}{f}}& f\&GreaterEqual;\stackrel{\&OverBar;}{f}\\ k_4& f<\stackrel{\&OverBar;}{f}\end{array}\right.---\left(17\right)$

式中：f_max是群体中最大的适应值，

是群体的平均适应值，f是个体的适应值，k₂，k₄是(0，1)间的常数，P_m是变异概率；

(5)将新老种群合并，进行适应度共享调整，

适值共享值调整的函数为：

$f_i^{\left(\mathrm{sh}\right)}=\frac{f_i}{m_i}=\frac{f_i}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\mathrm{sh}\left(d_{\mathrm{ij}}\right)}---\left(18\right)$

式中：N是种群个体数；m_i是个体i的小生境数；d_ij是个体i，j间的距离；sh(d_ij)是共享值，描述个体之间的相似度；

共享值的函数为：

$\mathrm{sh}\left(d_{\mathrm{ij}}\right)=\left\{\begin{array}{cc}1-{(d_{\mathrm{ij}}/{\mathrm{\&sigma;}}_{\mathrm{sh}})}^{\mathrm{\&alpha;}}& d_{\mathrm{ij}}<{\mathrm{\&sigma;}}_{\mathrm{sh}}\\ 0& d_{\mathrm{ij}}\&GreaterEqual;{\mathrm{\&sigma;}}_{\mathrm{sh}}\end{array}\right.---\left(19\right)$

式中：σ_sh是小生境半径；α是常数一般取1；

(6)对每一个个体，按照适值共享排挤策略，淘汰相似个体中适应度最差的，完成重插入，形成子代种群；

适值共享排挤的实现步骤为：

①把父代和子代放在一起形成新的种群，适应值共享调整群体的适应度；

②在父代个体中随机生成Cf个排挤因子组，每一组由随机生成的s个父代个体组成；

③在每一个排挤因子组中，寻找一个与子代个体最为相似的个体，这样一共可选出Cf个个体；

④将Cf+1个个体中适值最小的个体淘汰；

⑤重复步骤②～步骤④；

其中Cf的取值区间为[2，4]，而s的取值，至少为要搜索的峰个数的两倍；

(7)终止判断，如果满足终止条件，循环结束；如果不满足，将(6)的种群作为父代执行(2)～(6)。

10)执行下一层，重复步骤5)～9)，直至完成整个模型各向异性参数和界面的3次系数的优化。

本发明是采用“从浅到深剥层”的方式逐层计算，对于3D VSP情况，将其重排为不同方向的Walkaway VSP测线，同样适用。

本发明有如下特点：

模型为二维起伏层状各向异性VTI介质，反射界面采用三次函数表示更接近真实地层；

采用线性收敛的试射射线追踪计算效率高，简化了各向异性介质复杂的内在关系；

自适应交叉概率、变异概率避免了早熟、提高收敛速度；

适值共享排挤策略有利保持种群多样性，能有效处理复杂多峰问题；

利用Walkaway VSP初至旅行时、极化角建立反演目标函数，在遗传算法导向搜索指导下迭代试射射线追踪正演，非线性反演模型各向异性参数、界面形态。

附图说明

图1利用Walkaway VSP数据反演各向异性参数技术流程。

图2二维层状各向异性模型示意图。

图3VTI介质波前面、波射线、波矢量示意图。

图4VTI介质透反射关系示意图。

图5实例1模型参数及3个共深度点Walkaway VSP射线路径。

图6实列1的初至旅行时(上)、极化角(下)。

图7实例1初至旅行时反演模型初至与实测初至对比，从左到右依次为层2，层3，层4。

图8实例1极化角反演模型极化角与实测极化角对比。从左到右依次为层2，层3，层4。

图9实例1初至旅行时和极化角联合反演模型初至与实测初至对比。从左到右依次为层2，层3，层4。

图10实例1初至旅行时和极化角联合反演模型极化角与实测极化角对比。从左到右依次为层2，层3，层4。

图11实例2模型参数及4个共深度点Walkaway VSP射线路径。

图12实列2初至旅行时(上)、极化角(下)。

图13实列2反演值与真值对比。

图14野外实测零偏纵横波初至。

图15垂向纵横波速度及比值。

图16反演初始模型及观测系统。

图17初始模型初至与实测初至对比。

图18反演的各向异性参数及地层倾角。

图19反演的地层模型。

图20反演的模型初至与实测初至对比。

具体实施方式

本发明是利用变偏移距垂直地震剖面数据，采用“从浅到深剥层”的方式逐层反演各向异性参数优化界面形态；对于三维垂直地震剖面情况，重排为不同方向的变偏移距垂直地震剖测线，同样适用。

以下结合附图和3个实例进一步说明本发明的实施步骤和技术流程(图1)。

本发明采用以下技术步骤：

1)由变偏移距垂直地震剖面野外观测系统确定以井口为中心的模型范围，假定条件是各向异性模型为横向各向同性(VTI)介质；

2)用零偏垂直地震剖面(VSP)纵、横波初至旅行时计算垂向纵、横波速度，根据纵波速度的变化趋势划分地层并给出各层与井轨迹的交点，计算各层的层间垂向纵、横波速度；计算垂向横波速度时，若无零偏横波资料，则利用转换横波资料。

3)在变偏移距垂直地震剖面炮线方向的深度域过井地震剖面上，拾取标志层的层位，计算各层与井轨迹的交点；由距离最近的上下两个标志层，采用线性加权计算得到中间层界面的坐标(x_i，z_i)，由这些坐标拟合出界面的3次多项式函数；所述的线性加权函数是：

$z_i=\frac{(z^i-z_m^j)}{(z_m^{j+1}-z_m^j)}{z}_{m,j}+\frac{(z_m^{j+1}-z^i)}{(z_m^{j+1}-z_m^j)}{z}_{m,j+1}---\left(1\right)$

x_i＝x_m，j

式中：(x_i，z_i)是中间层i的坐标；i是地层号；j是标志层号；(xⁱ，zⁱ)是层i与井轨迹交点的坐标；

是标志层与井轨迹交点的纵坐标；

z_m，j、z_m，j+1是拾取的标志层纵坐标；

x_m，j是拾取的标志层的横坐标；

由中间层i的坐标(x_i，z_i)拟合出界面的3次多项式函数为：

z＝a₁·x³+a₂·x²+a₃·x+a₄    (2)

式中：a₁、a₂、a₃、a₄是界面的3次系数；

所建立的二维层状各向异性模型如图2；

4)由已知的表层速度资料建立表层模型，根据变偏移距垂直地震剖面野外观测情况给模型建立相应的观测系统；若无表层速度资料，则利用已知的静校正数据，建立不带表层结构的从某一基准面开始的模型，并根据变偏移距垂直地震剖面野外观测情况给模型建立相应的观测系统；

5)在参数的范围内，随机初始化待求层的各向异性参数和待求层顶界面的3次系数，形成待优化的各向异性速度模型；

6)根据待求层顶界面、底界面与井轨迹的交点，筛选位于该层之间的检波点，拾取这些检波点接收到得各个炮点的初至、计算极化角；

7)试射射线追踪正演步骤5)模型在步骤6)的检波点位置接收到的各个炮点初至、极化角；极化角是射线到达检波点时，射线方向与垂向所成的角度。试射射线追踪正演初至旅行时、极化角采用以下方法：

(1)给定初始相角；

(2)计算相速度、群速度、群角；

由相角计算相速度的公式为：

$V_P\left(\mathrm{\&theta;}\right)=V_{P0}\sqrt{1+\mathrm{\&epsiv;}{\sin }^2\mathrm{\&theta;}+\frac{1}{2}(1-\frac{V_{S0}^2}{V_{P0}^2})\left\{\right[1+\frac{4{\mathrm{\&delta;}}^{*}}{(1-V_{S0}^2/V_{P0}^2)}s{\mathrm{in}}^2\mathrm{\&theta;}{\cos }^2\mathrm{\&theta;}+\frac{4(1-V_{S0}^2/V_{P0}^2+\mathrm{\&epsiv;})\mathrm{\&epsiv;}}{(1-V_{S0}^2/V_{P0}^2)}{\sin }^4\mathrm{\&theta;}{]}^{1/2}-1\}}---\left(3\right)$

式中：θ是相角，ε、δ^*是Thomsen各向异性参数，V_P0、V_S0是纵、横波垂直各向同性面的相速度，V_P(θ)是相角θ对应的纵波相速度。

由相角计算群速度的公式为：

$V_g\left(\mathrm{\&phi;}\right)=\sqrt{V^2\left(\mathrm{\&theta;}\right)+{\left(\frac{\mathrm{dV}\left(\mathrm{\&theta;}\right)}{\mathrm{d\&theta;}}\right)}^2}---\left(4\right)$

式中：θ是相角，φ是群角，V(θ)是相速度，V_g(φ)是群速度。

由相速度、群速度计算群角的公式为：

$\mathrm{\&phi;}=\arccos \left(\frac{V_g\left(\mathrm{\&phi;}\right)}{V\left(\mathrm{\&theta;}\right)}\right)+\mathrm{\&theta;}---\left(5\right)$

式中：θ是相角，φ是群角，V(θ)是相速度，V_g(φ)是群速度。

(3)射线沿群角方向传播，计算射线与界面的交点；

(4)在界面交点处由斯奈尔定律做透射处理，计算透射相角；

所述的斯奈尔定律为：

$P=\frac{\sin \left({\mathrm{\&alpha;}}_{\mathrm{inc}}\right)}{V\left({\mathrm{\&theta;}}_{\mathrm{inc}}\right)}=\frac{\sin \left({\mathrm{\&alpha;}}_{\mathrm{tra}}\right)}{V\left({\mathrm{\&theta;}}_{\mathrm{tra}}\right)}---\left(6\right)$

式中：α_dip是界面倾角；θ_inc是入射相角，α_inc是入射角；θ_tra是透射相角，α_tra是透射角；

界面交点处的各个角度关系为：

α_inc＝|α_dip+θ_inc|

θ_tra＝sign(α_dip+θ_inc)α_tra-α_dip    (7)

式中：α_dip是界面倾角；θ_inc是入射相角，α_inc是入射角；θ_tra是透射相角，α_tra是透射角；

(5)重复(2)-(4)逐层计算直到射线到达井轨迹；

(6)计算射线与井轨迹的交点，计算交点与检波点的距离，若交点与检波点的距离小于给定精度，该射线即为给定精度条件下的准确射线，计算极化角；

若交点与检波点的距离大于给定精度，则由检波点附近的两条射线线性插值计算初始相角，重复(2)-(6)；

线性插值计算初始相角的公式为：

$\mathrm{\&alpha;}={\mathrm{\&alpha;}}_1+\frac{z-z_1}{z_2-z_1}({\mathrm{\&alpha;}}_2-{\mathrm{\&alpha;}}_1)---\left(8\right)$

式中：α₁、α₂是检波点附近的两条射线的初始相角；

(x₁，z₁)、(x₂，z₂)是α₁、α₂对应出射坐标(即射线与井轨迹的交点)；

(x，z)是检波点坐标，α是待求的初始相角。

(7)计算射线路径对应的旅行时；

由射线路径计算旅行时的公式为：

$T_{\mathrm{tra}}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{l_i}{V_{\mathrm{gi}}}---\left(9\right)$

式中，n是射线段数；

l_i是射线追踪到的第i个射线段；

V_gi是第i个射线段对应的群速度；

T_tra实现路径对应的旅行时；

VTI介质波前面、波射线、波矢量、群角、相角的关系如图3所示，透反射关系如图4所示；

8)采用P波初至旅行时或P波极化角或联合使用P波初至旅行时和P波极化角建立目标函数；采用P波初至旅行时建立的目标函数为：

$\mathrm{objv}=\sqrt{\frac{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{(T_{\mathrm{for}}^k-T_{\mathrm{med}}^k)}^2}{n}}---\left(10\right)$

式中：

n是实测的初至个数；

是该层检波点野外实测的第k个初至；

是该层检波点接收到的模型正演的第k个初至。

采用P波极化角建立的目标函数为：

$\mathrm{objv}=\sqrt{\frac{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&alpha;}}_{\mathrm{for}}^k-{\mathrm{\&alpha;}}_{\mathrm{med}}^k)}^2}{n}}---\left(11\right)$

式中：

n是实测的初至个数；

是该层检波点野外实测的第k个极化角；

是该层检波点接收到模型正演的第k个极化角。

步骤8)所述的联合使用P波初至旅行时和P波极化角建立的目标函数为：

$\mathrm{objv}=C_T\sqrt{\frac{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{(T_{\mathrm{for}}^k-T_{\mathrm{med}}^k)}^2}{n}}+C_{\mathrm{\&alpha;}}\sqrt{\frac{\underset{k=1}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&alpha;}}_{\mathrm{for}}^k-{\mathrm{\&alpha;}}_{\mathrm{med}}^k)}^2}{n}}---\left(12\right)$

式中：

n是实测的初至个数；

该层检波点野外实测的第k个初至；

该层检波点接收到模型正演的第k个初至；

该层检波点野外实测的第k个极化角；

该层检波点接收到模型正演的第k个极化角；

C_T、C_α是初至、极化角在目标函数中的权系数。

9)利用小生境遗传算法优化步骤8)的目标函数；

步骤8)目标函数最小值对应的步骤5)的待优化的各向异性速度模型就是最优模型，这就得到了待求层的层各向异性参数和该层顶界面的3次系数。

步骤9)中所述的小生境遗传算法的参数是：

编码：实数编码；

初始种群生成：在变量范围，随机生成；

选择算子：随机遍历采样选择；

交叉操作：均匀算术交叉，自适应交叉概率；

变异操作：实值变异，自适应变异概率；

群体间的迭代策略-适值共享排挤：将父代+子代形成新的种群，适值共享调整个体的适值，然后执行排挤操作直至种群个体数与父代相同；

终止条件：最大遗传代数+目标函数小于给定精度。

步骤9)中所述的小生境遗传算法的实现步骤为：

(1)随机产生初始化群体，对个体适应度进行评价，

个体适应度评价的函数为：

$\mathrm{FitnV}=2-\mathrm{sp}+2\&times;(\mathrm{sp}-1)\&times;\frac{\mathrm{Pos}-1}{\mathrm{Nind}-1}---\left(13\right)$

式中，sp是压差数；

Pos个体在种群中降序排序的位置；

Nind种群中个体的个数；

FitnV个体适应度；

(2)随机遍历采样选择配对父代个体，

随机遍历采样选择的概率为：

$F\left(x_i\right)=\frac{f\left(x_i\right)}{\underset{i=1}{\overset{\mathrm{Nind}}{\mathrm{\&Sigma;}}}f\left(x_i\right)}---\left(14\right)$

式中，f(x_i)是个体x_i的适应度；

Nind种群中个体的个数；

F(x_i)是个体x_i被选择的概率；

(3)按照自适应概率，均匀算术交叉选择配对的父代个体，

自适应交叉概率为：

$P_c=\left\{\begin{array}{cc}\frac{k_1(f_{\max }-f^{\&prime;})}{f_{\max }-\stackrel{\&OverBar;}{f}}& f^{\&prime;}\&GreaterEqual;\stackrel{\&OverBar;}{f}\\ k_3& f^{\&prime;}<\stackrel{\&OverBar;}{f}\end{array}\right.---\left(15\right)$

式中：f_max是群体中最大的适应值，

是群体的平均适应值，f′是交叉的两个个体中较大的适应值，k₁，k₃是(0，1)间的常数，P_c是交叉概率；

均匀算术交叉函数为：

$X_A^{t+1}=\mathrm{\&alpha;}{X}_B^t+(1-\mathrm{\&alpha;})X_A^t$ (16)

$X_B^{t+1}=\mathrm{\&alpha;}{X}_A^t+(1-\mathrm{\&alpha;})X_B^t$

式中：

是待交叉的父代个体；

α是[0，1]之间的一个常数；

是交叉后的子代个体；

(4)按照自适应概率，实值变异交叉过的个体，产生新的种群，

自适应变异概率为：

$P_m=\left\{\begin{array}{cc}\frac{k_2(f_{\max }-f)}{f_{\max }-\stackrel{\&OverBar;}{f}}& f\&GreaterEqual;\stackrel{\&OverBar;}{f}\\ k_4& f<\stackrel{\&OverBar;}{f}\end{array}\right.---\left(17\right)$

式中：f_max是群体中最大的适应值，是群体的平均适应值，f是个体的适应值，k₂，k₄是(0，1)间的常数，P_m是变异概率；

(5)将新老种群合并，进行适应度共享调整，

适值共享值调整的函数为：

$f_i^{\left(\mathrm{sh}\right)}=\frac{f_i}{m_i}=\frac{f_i}{\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}\mathrm{sh}\left(d_{\mathrm{ij}}\right)}---\left(18\right)$

式中：N是种群个体数；m_i是个体i的小生境数；d_ij是个体i，j间的距离；sh(d_ij)是共享值，描述个体之间的相似度；

共享值的函数为：

$\mathrm{sh}\left(d_{\mathrm{ij}}\right)=\left\{\begin{array}{cc}1-{(d_{\mathrm{ij}}/{\mathrm{\&sigma;}}_{\mathrm{sh}})}^{\mathrm{\&alpha;}}& d_{\mathrm{ij}}<{\mathrm{\&sigma;}}_{\mathrm{sh}}\\ 0& d_{\mathrm{ij}}\&GreaterEqual;{\mathrm{\&sigma;}}_{\mathrm{sh}}\end{array}\right.---\left(19\right)$

式中：σ_sh是小生境半径；α是常数一般取1；

(6)对每一个个体，按照适值共享排挤策略，淘汰相似个体中适应度最差的，完成重插入，形成子代种群；

适值共享排挤的实现步骤为：

①把父代和子代放在一起形成新的种群，适应值共享调整群体的适应度；

②在父代个体中随机生成Cf个排挤因子组，每一组由随机生成的s个父代个体组成；

③在每一个排挤因子组中，寻找一个与子代个体最为相似的个体，这样一共可选出Cf个个体；

④将Cf+1个个体中适值最小的个体淘汰；

⑤重复步骤②～步骤④；

其中Cf的取值区间为[2，4]，而s的取值，至少为要搜索的峰个数的两倍；

(7)终止判断，如果满足终止条件，循环结束；如果不满足，将(6)的种群作为父代执行(2)～(6)。

10)执行下一层，重复步骤5)～9)，直至完成整个模型各向异性参数和界面的3次系数的优化。

以下进一步通过实例说明本发明。

实例1理论模型图5，大小4000m×2000m，井位于X＝2000m，起始炮点X＝80m，炮点间距80m，49炮采集，检波点位于Z＝800m、1200m、1600m，初至旅行时及极化角见图6。模型包括4个反射界面，1个水平的、1个单斜的、2个起伏的。反演时假设已知垂向纵横波速度和反射界面，从浅层到深层依次计算，求2，3，4层的各层各向异性参数ε、δ^*。

(1)采用初至旅行时建立目标函数，第二层各向异性参数反演相对误差为-0.0003、-0.0031，旅行时标准差std＝0.006ms；第三层各向异性参数反演相对误差为-0.0028、0.0097，旅行时标准差std＝0.032ms；第四层各向异性反演参数相对误差为-0.0032、0.0015，旅行时标准差std＝0.0212ms。反演结果模型初至与实测初至对比，见图7。

(2)采用极化角建立目标函数，第二层各向异性参数反演相对误差为0.0000、0.0001，极化角标准差std＝7.9e-7rand；第三层各向异性参数反演相对误差为0.0001、-0.0001，极化角标准差std＝2.72e-6rand；第四层各向异性参数反演相对误差为-0.0002、0.0003，极化角标准差std＝1.18e-4rand。反演结果模型极化角与实测极化角对比，见图8。

(3)联合使用初至旅行和极化角建立目标函数，第二层各向异性参数反演相对误差为-0.0002、0.0024，旅行时标准差std＝0.0038ms，极化角标准差std＝3.3e-5rand；第三层各向异性参数反演相对误差为0.0002、-0.0019，旅行时标准差std＝0.0027ms，极化角标准差std＝5.5e-5rand；第四层各向异性参数反演相对误差为-0.0002、0.0001，旅行时标准差std＝0.0022ms，极化角标准差std＝1.3e-5rand。反演结果模型初至与实测初至对比，见图9；反演结果模型极化角与实测极化角对比，见图10。

实例2理论模型图11，模型大小4000m×3000m，井位于X＝2000m，起始炮点X＝80m，炮点间距80m，49炮采集，检波点位于Z＝700m、1200m、1900m、2400m，初至旅行时及极化角见图12。模型包括5个反射界面，1个水平的、1个单斜的、3个起伏的。反演时假设已知垂向纵横波速度，从浅层到深层依次计算，求层2、3、4、5的各向异性参数和反射界面1、2、3、4的形态参数。联合使用初至旅行时和极化角建立目标函数，反演结果误差分析见表1，反演值与真值对比见图13。

表1实例2反演结果误差分析

实例3野外实测WVSP资料，43炮采集，炮间距180m，16级采集，级间距10m，共观测160个深度点。在该井同时还进行了零偏纵波、零偏横波采集；根据零偏纵波初至、横波初至及其变化，见图14，分为22层，建立速度深度模型，同时得到各层的纵横波套层速度及纵横波速度比，见图15。各炮点以井口为参考变换到模型坐标，各层位、接收点均以井口为参考点从实际坐标映射到模型坐标。

从收集的资料来看，该WVSP观测范围为地层基本上可用斜层表示。建立平层初始速度深度模型，选取12个深度点用于各向异性参数反演和反射界面调整，见图16；初始模型初至与实测初至对比，见图17，可以发现模型地层的倾向为左低右高。表层当做一套层处理，WVSP最浅的观测点以上地层视为等效单斜层处理(即地层倾角相同，各向异性系数相同)。反演的各向异性参数、地层倾角，见图18；反演得到地层模型，见图19，从图20可以看到模型初至与实测初至吻合的很好。

Claims (11)

1.一种利用变偏移距垂直地震剖面数据反演各向异性参数的方法，特点是采用以下技术步骤：

1)由变偏移距垂直地震剖面野外观测系统确定以井口为中心的模型范围，假定条件是各向异性模型为横向各向同性(VTI)介质；

2)用零偏垂直地震剖面(VSP)纵、横波初至旅行时计算垂向纵、横波速度，根据纵波速度的变化趋势划分地层并给出各层与井轨迹的交点，计算各层的层间垂向纵、横波速度；

3)在变偏移距垂直地震剖面炮线方向的深度域过井地震剖面上，拾取标志层的层位，计算各层与井轨迹的交点；由距离最近的上下两个标志层，采用线性加权计算得到中间层界面的坐标(x_i，z_i)，由这些坐标拟合出界面的3次多项式函数；

4)由已知的表层速度资料建立表层模型，根据变偏移距垂直地震剖面野外观测情况给模型建立相应的观测系统；

5)在参数的范围内，随机初始化待求层的各向异性参数和待求层顶界面的3次系数，形成待优化的各向异性速度模型；

6)根据待求层顶界面、底界面与井轨迹的交点，筛选位于该层之间的检波点，拾取这些检波点接收到得各个炮点的初至、计算极化角；

7)试射射线追踪正演步骤5)模型在步骤6)的检波点位置接收到的各个炮点初至、极化角；

8)采用P波初至旅行时或P波极化角或联合使用P波初至旅行时和P波极化角建立目标函数；

9)利用小生境遗传算法优化步骤8)的目标函数；

10)执行下一层，重复步骤5)～9)，直至完成整个模型各向异性参数和界面的3次系数的优化。

2.根据权利要求1所述的方法，特点是步骤2)计算垂向横波速度时，若无零偏横波资料，则利用转换横波资料。

3.根据权利要求1所述的方法，特点是步骤3)所述的线性加权函数是：

x_i＝x_m，j

式中：(x_i，z_i)是中间层i的坐标；i是地层号；j是标志层号；(xⁱ，zⁱ)是层i与井轨迹交点的坐标；

是标志层与井轨迹交点的纵坐标；

z_m，j、z_m，j+1是拾取的标志层纵坐标；

x_m，j是拾取的标志层的横坐标；

由中间层i的坐标(x_i，z_i)拟合出界面的3次多项式函数为：

z＝a₁·x³+a₂·x²+a₃·x+a₄    (2)

式中：a₁、a₂、a₃、a₄是界面的3次系数。

4.根据权利要求1所述的方法，特点是步骤4)若无表层速度资料，则利用已知的静校正数据，建立不带表层结构的从某一基准面开始的模型，并根据变偏移距垂直地震剖面野外观测情况给模型建立相应的观测系统。

5.根据权利要求1所述的方法，特点是步骤7)所述的极化角是射线到达检波点时，射线方向与垂向所成的角度。

6.根据权利要求1所述的方法，特点是步骤7)所述的试射射线追踪正演初至旅行时、极化角采用以下方法：

(1)给定初始相角；

(2)计算相速度、群速度、群角；

由相角计算相速度的公式为：

式中：θ是相角，ε、δ^*是Thomsen各向异性参数，V_P0、V_S0是纵、横波垂直各向同性面的相速度，V_P(θ)是相角θ对应的纵波相速度。

由相角计算群速度的公式为：

式中：θ是相角，φ是群角，V(θ)是相速度，V_g(φ)是群速度。

由相速度、群速度计算群角的公式为： 

式中：θ是相角，φ是群角，V(θ)是相速度，V_g(φ)是群速度。

(3)射线沿群角方向传播，计算射线与界面的交点；

(4)在界面交点处由斯奈尔定律做透射处理，计算透射相角；

所述的斯奈尔定律为：

式中：α_dip是界面倾角；θ_inc是入射相角，α_inc是入射角；θ_tra是透射相角，α_tra是透射角；

界面交点处的各个角度关系为：

α_inc＝|α_dip+θ_inc|

θ_tra＝sign(α_dip+θ_inc)α_tra-α_dip    (7)

式中：α_dip是界面倾角；θ_inc是入射相角，α_inc是入射角；θ_tra是透射相角，α_tra是透射角；

(5)重复(2)-(4)逐层计算直到射线到达井轨迹；

(6)计算射线与井轨迹的交点，计算交点与检波点的距离，若交点与检波点的距离小于给定精度，该射线即为给定精度条件下的准确射线，计算极化角；

若交点与检波点的距离大于给定精度，则由检波点附近的两条射线线性插值计算初始相角，重复(2)-(6)；

线性插值计算初始相角的公式为：

式中：α₁、α₂是检波点附近的两条射线的初始相角；

(x₁，z₁)、(x₂，z₂)是α₁、α₂对应出射坐标(即射线与井轨迹的交点)；

(x，z)是检波点坐标，α是待求的初始相角。

(7)计算射线路径对应的旅行时；

由射线路径计算旅行时的公式为： 

式中，n是射线段数；

l_i是射线追踪到的第i个射线段；

V_gi是第i个射线段对应的群速度；

T_tra实现路径对应的旅行时。

7.根据权利要求1所述的方法，特点是步骤8)所述的采用P波初至旅行时建立的目标函数为：

式中：

n是实测的初至个数；

是该层检波点野外实测的第k个初至；

是该层检波点接收到的模型正演的第k个初至。

8.根据权利要求1所述的方法，特点是步骤8)所述的采用P波极化角建立的目标函数为：

式中：

n是实测的初至个数；

是该层检波点野外实测的第k个极化角；

是该层检波点接收到模型正演的第k个极化角。

9.根据权利要求1所述的方法，特点是步骤8)所述的联合使用P波初至旅行时和P波极化角建立的目标函数为：

式中：

n是实测的初至个数；

该层检波点野外实测的第k个初至；

该层检波点接收到模型正演的第k个初至；

该层检波点野外实测的第k个极化角；

该层检波点接收到模型正演的第k个极化角；

C_T、C_α是初至、极化角在目标函数中的权系数。

10.根据权利要求1所述的方法，特点是步骤8)目标函数最小值对应的步骤5)的待优化的各向异性速度模型就是最优模型，这就得到了待求层的层各向异性参数和该层顶界面的3次系数。

11.根据权利要求1所述的方法，特点是步骤9)中所述的小生境遗传算法的参数是：

编码：实数编码；

初始种群生成：在变量范围，随机生成；

选择算子：随机遍历采样选择；

交叉操作：均匀算术交叉，自适应交叉概率；

变异操作：实值变异，自适应变异概率；

群体间的迭代策略-适值共享排挤：将父代+子代形成新的种群，适值共享调整个体的适值，然后执行排挤操作直至种群个体数与父代相同；

终止条件：最大遗传代数+目标函数小于给定精度。

步骤9)中所述的小生境遗传算法的实现步骤为：

(1)随机产生初始化群体，对个体适应度进行评价，

个体适应度评价的函数为：

式中，sp是压差数；

Pos个体在种群中降序排序的位置；

Nind种群中个体的个数；

FitnV个体适应度；

(2)随机遍历采样选择配对父代个体， 

随机遍历采样选择的概率为：

式中，f(x_i)是个体x_i的适应度；

Nind种群中个体的个数；

F(x_i)是个体x_i被选择的概率；

(3)按照自适应概率，均匀算术交叉选择配对的父代个体，

自适应交叉概率为：

式中：f_max是群体中最大的适应值， 

是群体的平均适应值，f′是交叉的两个个体中较大的适应值，k₁，k₃是(0，1)间的常数，P_c是交叉概率；

均匀算术交叉函数为：

(16)

式中： 

是待交叉的父代个体；

α是[0，1]之间的一个常数；

是交叉后的子代个体；

(4)按照自适应概率，实值变异交叉过的个体，产生新的种群，

自适应变异概率为：

式中：f_max是群体中最大的适应值， 

是群体的平均适应值，f是个体的适应值，k₂，k₄是(0，1)间的常数，P_m是变异概率； 

(5)将新老种群合并，进行适应度共享调整，

适值共享值调整的函数为：

式中：N是种群个体数；m_i是个体i的小生境数；d_ij是个体i，j间的距离；sh(d_ij)是共享值，描述个体之间的相似度；

共享值的函数为：

式中：σ_sh是小生境半径；α是常数一般取1；

(6)对每一个个体，按照适值共享排挤策略，淘汰相似个体中适应度最差的，完成重插入，形成子代种群；

适值共享排挤的实现步骤为：

①把父代和子代放在一起形成新的种群，适应值共享调整群体的适应度；

②在父代个体中随机生成Cf个排挤因子组，每一组由随机生成的s个父代个体组成；

③在每一个排挤因子组中，寻找一个与子代个体最为相似的个体，这样一共可选出Cf个个体；

④将Cf+1个个体中适值最小的个体淘汰；

⑤重复步骤②～步骤④；

其中Cf的取值区间为[2，4]，而s的取值，至少为要搜索的峰个数的两倍；

(7)终止判断，如果满足终止条件，循环结束；如果不满足，将(6)的种群作为父代执行(2)～(6)。 

Images