本发明公开了一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，该混合反演策略将线性和非线性反演相结合，利用线性反演结果作为约束，生成较小的搜索范围和初始粒子，开展非线性反演以同时获取五参数的反演结果。该方法克服常规线性反演易落入局部极值的缺点，提高非线性反演的计算效率和准确性。采用反射系数精确式可克服近似表达式小角度、弱反射和弱各向异性的限制，提高反演的适用性。测井模型数据测试验证了方法的有效性，这为后续VTI介质各向异性研究奠定了基础。
1.一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，其特征在于，包括以下步骤：
步骤一：获取PP波叠前角度道集，并从中提取地震子波，并建立用于反演的子波矩阵；
步骤二：获取工区的井孔数据和三维层位解释结果，建立反演所用的初始模型；
步骤三：预设线性反演循环迭代的最大迭代次数，开始线性反演循环迭代，利用初始模型和子波矩阵获取合成地震记录；
步骤四：利用初始模型计算反射系数精确式正演算子的偏导数；
步骤五：建立目标函数，基于正演算子的偏导数和合成地震记录，计算目标函数的雅可比矩阵和伪海森矩阵；
步骤六：利用目标函数的雅可比矩阵和伪海森矩阵计算目标参数向量的更新方向，利用迭代格式更新目标参数向量；
步骤七：若线性反演循环迭代次数未达到预设的线性反演循环迭代的最大迭代次数，则开始下一次线性反演迭代，重复步骤三至六；若达到，则结束线性反演循环迭代过程，输出线性反演结果，继续执行步骤八；
步骤八：利用线性反演结果，生成搜索范围和初始粒子；
步骤八中，针对目标参数向量中某一深度样点，用m 表示线性反演结果，则基于线性反min max
其中，m 为搜索范围的最小值，m 为搜索范围的最大值，m 表示井口处抽取到的真实参考值，b为尺度因子；
设定N个初始粒子，粒子群x表示为：其中，j的范围是1到N，xj表示第j个粒子，也是第j个可能的模型参数；并将线性反演结Li果m 赋值给N个粒子中的任意一个，最终形成初始粒子群；
步骤九：预设非线性反演迭代的最大迭代次数，开始非线性反演循环迭代，将目标函数应用于所有初始粒子计算适应度值，并从所有初始粒子的适应度值中挑选出个体最优位置步骤十：将目标函数应用于所有粒子来计算替换适应度值，并计算不同粒子的接收概gbest
步骤十一：根据个体最优位置和群体最优位置来计算更新速度，进而采用迭代格式更新粒子；
步骤十二：若迭代次数未达到预设的非线性反演迭代的最大迭代次数，则开始下一次非线性反演循环迭代，重复步骤九至十一；若达到最大迭代次数，则结束非线性反演，输出所有粒子，选取所有粒子中的最好结果作为最终的非线性反演结果。
2.根据权利要求1所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，其特征在于，步骤二中，井孔数据包括纵波速度α、横波速度β、密度ρ，各向异性参数ε和δ曲线，与四个刚度参数c11、c13、c33和c55之间的关系为：
c11＝(2ε+1)αρ,c13=ρη‑βρ其中，η为中间变量；
3.根据权利要求2所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，其特征在于，步骤三中，获取合成地震记录具体如下：利用VTI介质反射系数精确方程计算反射系数r：A·r＝b
其中，A和b均为中间过程矩阵，具有如下表达形式：其中， 为上层介质纵波相关的中间变量， 为上层介质横波相关的中间变量， 为下层介质纵波相关的中间变量， 为下层介质横波相关的中间变量；
为上层介质的刚度参数， 为下层介质的刚度参数； 分别为上、下层纵波相关的垂直慢度， 分别为上、下层横波相关的垂直慢度，p为水平慢度；
p这些参数均由刚度参数c11、c13、c33、c55和密度ρ计算获得，上标T表示转置运算符；
基于褶积理论，合成地震记录dsyn利用子波矩阵W和反射系数r计算求得：dsyn＝G(m)＝W·r
m＝(α1，…αi，…αM，β1，…βi，…βM，ρ1，…ρi，…ρM，ε1，…εi，…εM，δ1，…δi，…δM)αi、βi、ρi、εi和δi分别表示为第i个时间采样点所对应的纵波速度、横波速度、密度、两个各向异性参数，i＝1,....,M；M为模型时间深度最大样点个数；G(m)为反射系数精确式正演算子，即合成地震记录。
4.根据权利要求3所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，其特征在于，步骤四中，计算反射系数精确式正演算子的偏导数具体如下：反射系数精确式正演算子G(m)对目标参数向量m的偏导数为：其中， 为反射系数对目标参数向量的偏导数， 包括对反射界面上层介质模型参U L
数m和下层介质模型参数m的偏导数 具体的计算公式为：U
式中， 分别表示中间变量矩阵A对反射界面上层介质模型参数m 和下层介L
质模型参数m 的偏导数， 分别表示中间变量矩阵b对反射界面上层介质模型参U L
5.根据权利要求4所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，其特征在于，步骤五中，建立目标函数S(m)，表示为：其中，d为实际叠前地震道集数据，λ为正则化参数，Cm为m计算得到的协方差矩阵，u为m的期望；
计算目标函数S(m)的雅可比矩阵J(m)，即一阶偏导数矩阵：其中， 为G(m)对m的偏导数；
针对目标函数S(m)的伪海森矩阵H(m)，采用了内部小循环迭代形式来计算获得：对于线性反演循环迭代中第k次迭代、内部小循环首次迭代的伪海森矩阵Hk，1(m)，由下式计算得到：对于线性反演循环迭代中第k次迭代、内部小循环第i次迭代的为海森矩阵Hk，i(m)，由下式计算获得：
Hk，i(m) ＝(Qk，i‑1) Hk，i‑1(m) (Qk，i‑1)+Kk，i‑1sk，i‑1(sk，i‑1)其中，Qk，i‑1、Kk，i‑1、sk，i‑1为线性反演循环中第k次迭代、内部小循环第i‑1次迭代的中间过程矩阵，Hk，i‑1(m)为线性反演循环中第k次迭代内部小循环第i‑1次迭代的伪海森矩阵。
6.根据权利要求5所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，其特征在于，步骤六中，利用目标函数的雅可比矩阵和伪海森矩阵计算目标参数向量更新方向，利用迭代格式更新目标参数向量，具体形式如下：线性反演循环迭代的第k次迭代时，目标参数向量的更新方向‑gk(m)由下式计算得到：‑gk(m)＝‑Hk(m)·Jk(m)其中，Hk(m)和Jk(m)分别代表线性反演循环中第k次迭代所对应的伪海森矩阵和雅可比矩阵；
利用如下的更新迭代格式来更新mk+1：mk+1＝mk+vk+1，
其中，mk+1、mk分别表示线性反演循环中第k+1次、第k次迭代的目标参数向量，vk、vk+1分别表示线性反演循环中第k次、第k+1次迭代的中间过程矩阵，是动量系数，a为更新步长。
7.根据权利要求6所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，其特征在于，步骤九具体如下：
将第j个粒子xj代入目标函数计算得到S(xj)，S(xj)也被称为第j个粒子的适用度；计算出所有粒子的适用度后，选择适用度最小的粒子赋值为第j个粒子的个体最优位置当l＝1,
其中，l为非线性反演循环迭代次数，xj,l表示第l次迭代中的第j个粒子，S(xj,l)、S(xj,l‑1)分别为第l次、第l‑1次迭代中的第j个粒子的适应度。
8.根据权利要求7所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，其特征在于，步骤十中，将目标函数应用于所有粒子来计算替换适应度值，并计算不同粒子gbest的接收概率，利用轮盘赌获取群体最优位置x 包括：gbest
计算替换适用度值Zr(xj,x ,Tl)：gbest gbest gbest其中，xj表示第j个粒子，x 为群体最优位置，S(xj)和S(x )分别表示将xj和xgbest
代入目标函数计算得到的适用度， 为xj和x 适用度的差；l为非线性反演循环中的迭代次数，Tl为第l次迭代所对应的温度值，τ表示为一个常数；
计算好所有粒子的替换适用度值后，计算不同粒子的接收概率Pg:其中，∑为求和运算符；
在得到所有粒子的接收概率后，采用轮盘赌的方式获得群体最优位置x 。
9.根据权利要求8所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，其特征在于，步骤十一中，根据个体最优位置和群体最优位置来计算更新速度，进而采用迭代格式更新粒子包括：gbest
基于已获得个体最优位置 和群体最优位置x ，计算更新速度：其中，vj,l+1、vj,l表示第j个粒子在第l+1次、第l次迭代中的更新速度，rand1(·)和rand2(·)表示两个介于[0,1]的随机数；κ2和κ3为两个控制更新步长的常数；xj,l表示第l次迭代的第j个粒子； 表示第l次迭代的第j个粒子的个体最优位置；
基于更新速度，更新不同的粒子：xj,l+1＝xj,l+κ1·vj,l+1其中，κ1为用以控制更新的步长的常数；xj,l+1表示第l+1次迭代的第j个粒子。
一种基于反射系数精确式的页岩VTI储层叠前混合反演方法
技术领域
[0001] 本发明涉及非常规各向异性储层地震勘探技术领域，特别是一种基于反射系数精确式的页岩VTI储层叠前混合反演方法。
背景技术
[0002] 叠前地震反演通过挖掘地震数据中的AVO/AVA(振幅随偏移距/入射角度变化)响应来获取地下介质的弹性信息。目前，叠前反演已是地震勘探的重要技术。以往，针对常规砂泥岩储层，均采用基于各向同性假设的反演方法，并在过去几十年里取得了巨大的成功。但随着油气需求不断增加，勘探的重点逐渐由常规储层向非常规油气转变。其中，页岩油气已经成为了目前勘探的重点与难点。页岩岩石表现为较强的VTI各向异性特征。VTI介质，即具有垂直对称轴的横向各向同性介质(Transverse isotropy with vertical axis of symmetry,VTI)是地下介质中存在最普遍的各向异性。以往基于各向同性假设的地震技术已不再适用于页岩储层的勘探。因此，为提高页岩储层的勘探精度，需开展适用于VTI介质的反演方法研究，用以有效挖掘页岩储层的弹性及各向异性信息。
[0003] 目前，针对VTI介质的正演研究较为完善，但是针对VTI介质的叠前反演方法并不完善，主要是选用反射系数近似式为正演算子。侯栋甲等(2014)将反射系数近似式引入VTI五参数直接反演中，并采用贝叶斯方法建立目标函数。但由于各向异性介质反演需要同时提取五个目标参数，且各向异性参数的地震敏感性较低，这导致VTI反演的不适定性更强，很难同时获取较好的五参数结果，尤其是各向异性参数的估算一直不理想。因此，Zhang等(2019)基于近似式提出了分步反演策略，通过数学转换将五参数转换为其他三项目标参数，然后通过间接反演方法换算得到各向异性参数，但是该方法引入了额外假设，且仍采用精度较低的近似式。VTI介质反射系数精确表达式虽然已经被给出，虽然其复杂的数学表达式，导致偏导数求取较难，因此目前较少被应用于叠前反演中。
[0004] 常规线性反演易落入局部极值，而非线性反演存在搜索范围大计算效率低的问题，近似式反演远角度精度低、弱各向异性限制等问题。
发明内容
[0005] 本发明所要解决的技术问题是克服现有技术的不足而提供一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，针对目前页岩等VTI(具有垂直同相轴的横向各向同性)介质的反演现存的难点问题，提出了一种将线性和非线性优化算法相结合的混合反演策略，并将VTI介质的反射系数精确式引入作为正演算子，用以提高页岩等VTI介质的各向异性参数的反演精度。
[0006] 本发明为解决上述技术问题采用以下技术方案：
[0007] 根据本发明提出的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法，包括以下步骤：
[0008] 步骤一：获取PP波叠前角度道集，并从中提取地震子波，并建立用于反演的子波矩阵；
[0009] 步骤二：获取工区的井孔数据和三维层位解释结果，建立反演所用的初始模型；
[0010] 步骤三：预设线性反演循环迭代的最大迭代次数，开始线性反演循环迭代，利用初始模型和子波矩阵获取合成地震记录；
[0011] 步骤四：利用初始模型计算反射系数精确式正演算子的偏导数；
[0012] 步骤五：建立目标函数，基于正演算子的偏导数和合成地震记录，计算目标函数的雅可比矩阵和伪海森矩阵；
[0013] 步骤六：利用目标函数的雅可比矩阵和伪海森矩阵计算目标参数向量的更新方向，利用迭代格式更新目标参数向量；
[0014] 步骤七：若线性反演循环迭代次数未达到预设的线性反演循环迭代的最大迭代次数，则开始下一次线性反演迭代，重复步骤三至六；若达到，则结束线性反演循环迭代过程，输出线性反演结果，继续执行步骤八；
[0015] 步骤八：利用线性反演结果，生成搜索范围和初始粒子；
[0016] 步骤九：预设非线性反演迭代的最大迭代次数，开始非线性反演循环迭代，将目标函数应用于所有初始粒子计算适应度值，并从所有初始粒子的适应度值中挑选出个体最优位置
[0017] 步骤十：将目标函数应用于所有粒子来计算替换适应度值，并计算不同粒子的接gbest收概率，利用轮盘赌获取群体最优位置x ；
[0018] 步骤十一：根据个体最优位置和群体最优位置来计算更新速度，进而采用迭代格式更新粒子；
[0019] 步骤十二：若迭代次数未达到预设的非线性反演迭代的最大迭代次数，则开始下一次非线性反演循环迭代，重复步骤九至十一；若达到最大迭代次数，则结束非线性反演，输出所有粒子，选取所有粒子中的最好结果作为最终的非线性反演结果。
[0020] 作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案，步骤二中，井孔数据包括纵波速度α、横波速度β、密度ρ，各向异性参数ε和δ曲线，与四个刚度参数c11、c13、c33和c55之间的关系为：
[0021] c33＝α2ρ,c55＝β2ρ
[0022] c11＝(2ε+1)α2ρ,c13＝ρη‑β2ρ
[0023] 其中，η为中间变量；
[0024]
[0025] 作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案，步骤三中，获取合成地震记录具体如下：
[0026] 利用VTI介质反射系数精确方程计算反射系数r：
[0027] A·r＝b
[0028] 其中，A和b均为中间过程矩阵，具有如下表达形式：
[0029]
[0030]
[0031] 其中， 为上层介质纵波相关的中间变量， 为上层介质横波相关的中间变量， 为下层介质纵波相关的中间变量， 为下层介质横波相关的中间变
量； 为上层介质的刚度参数， 为下层介质的刚度参数； 分
别为上、下层纵波相关的垂直慢度， 分别为上、下层横波相关的垂直慢度，p为水平慢度； p这些参数均由刚度参数c11、
c13、c33、c55和密度ρ计算获得，上标T表示转置运算符；
[0032] 基于褶积理论，合成地震记录dsyn利用子波矩阵W和反射系数r计算求得：
[0033] dsyn＝G(m)＝W·r
[0034] 其中，m为设定的目标参数向量：
[0035] m＝(α1,...,αM,β1,...,βM,ρ1,...,ρM,ε1,...,εM,δ1,...,δM)T[0036] αi、βi、ρi、εi和δi分别表示为第i个时间采样点所对应的纵波速度、横波速度、密度、两个各向异性参数，i＝1,....,M；M为模型时间深度最大样点个数，G(m)为反射系数精确式正演算子，即合成地震记录。
[0037] 作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案，步骤四中，计算反射系数精确式正演算子的偏导数具体如下：
[0038] 反射系数精确式正演算子G(m)对目标参数向量m的偏导数为：
[0039]
[0040] 其中， 为反射系数对目标参数向量的偏导数， 包括对反射界面上层介质模U L型参数m和下层介质模型参数m的偏导数 具体的计算公式为：
[0041]
[0042]
[0043] 式中， 分别表示中间变量矩阵A对反射界面上层介质模型参数mU和下L
层介质模型参数m的偏导数， 分别表示中间变量矩阵b对反射界面上层介质模
U L
型参数m和下层介质模型参数m的偏导数。
[0044] 作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案，步骤五中，
[0045] 建立目标函数S(m)，表示为：
[0046] S(m)＝[d‑G(m)]T[d‑G(m)]+λ(m‑u)TC‑m1(m‑u)
[0047] 其中，d为实际叠前地震道集数据，λ为正则化参数，Cm为m计算得到的协方差矩阵，u为m的期望；
[0048] 计算目标函数S(m)的雅可比矩阵J(m)，即一阶偏导数矩阵：
[0049]
[0050] 其中， 为G(m)对m的偏导数；
[0051] 针对目标函数S(m)的伪海森矩阵H(m)，采用了内部小循环迭代形式来计算获得：
[0052] 对于线性反演循环迭代中第k次迭代、内部小循环首次迭代的伪海森矩阵H k,1(m)，由下式计算得到：
[0053]
[0054] 对于线性反演循环迭代中第k次迭代、内部小循环第i次迭代的为海森矩阵Hk,i(m)，由下式计算获得：
[0055] Hk,i(m)‑1＝(Qk,i‑1)THk,i‑1(m)‑1(Qk,i‑1)+Kk,i‑1sk,i‑1(sk,i‑1)T[0056] 其中，Qk,i‑1、Kk,i‑1、sk,i‑1为线性反演循环中第k次迭代、内部小循环第i‑1次迭代的中间过程矩阵，Hk,i‑1(m)为线性反演循环中第k次迭代内部小循环第i‑1次迭代的伪海森矩阵。
[0057] 作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案，步骤六中，利用目标函数的雅可比矩阵和伪海森矩阵计算目标参数向量更新方向，利用迭代格式更新目标参数向量，具体形式如下：
[0058] 线性反演循环迭代的第k次迭代时，目标参数向量的更新方向‑gk(m)由下式计算得到：
[0059] ‑gk(m)＝‑Hk(m)·Jk(m)
[0060] 其中，Hk(m)和Jk(m)分别代表线性反演循环中第k次迭代所对应的伪海森矩阵和雅可比矩阵；
[0061] 利用如下的更新迭代格式来更新mk+1：
[0062]
[0063] 其中，mk+1、mk分别表示线性反演循环中第k+1次、第k次迭代的目标参数向量，vk、vk+1分别表示线性反演循环中第k次、第k+1次迭代的中间过程矩阵，是动量系数，a为更新步长。
[0064] 作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进Li一步优化方案，步骤八中，针对目标参数向量中某一深度样点，用m 表示线性反演结果，则min max
基于线性反演结果生成搜索范围[m m ]:
Re Li Li
[0065] 当|m ‑m |＜b·mRe Li Li
[0066] 当|m ‑m |≥b·m
[0067] 其中，mmin为搜索范围的最小值，mmax为搜索范围的最大值，mRe表示井口处抽取到的真实参考值，b为尺度因子；
[0068] 设定N个初始粒子，粒子群x表示为：
[0069]
[0070] 其中，j的范围是1到N，xj表示第j个粒子，也是第j个可能的模型参数；并将线性反Li演结果m 赋值给N个粒子中的任意一个，最终形成初始粒子群。
[0071] 作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案，步骤九具体如下：
[0072] 将第j个粒子xj代入目标函数计算得到S(xj)，S(xj)也被称为第j个粒子的适用度；计算出所有粒子的适用度后，选择适用度最小的粒子赋值为第j个粒子的个体最优位置[0073] 当l＝1,
[0074] 当l>1,
[0075] 其中，l为非线性反演循环迭代次数，xj,l表示第l次迭代中的第j个粒子，S(xj,l)、S(xj,l‑1)分别为第l次、第l‑1次迭代中的第j个粒子的适应度。
[0076] 作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案，步骤十中，将目标函数应用于所有粒子来计算替换适应度值，并计算不同粒gbest子的接收概率，利用轮盘赌获取群体最优位置x 包括：
[0077] 计算替换适用度值Zr(xj,xgbest,Tl)：
[0078]
[0079]
[0080] 其中，xj表示第j个粒子，xgbest为群体最优位置，S(xj)和S(xgbest)分别表示将xj和gbest gbestx 代入目标函数计算得到的适用度， 为xj和x 适用度的差；l为非线性反演循环中的迭代次数，Tl为第l次迭代所对应的温度值，τ表示为一个常数；
[0081] 计算好所有粒子的替换适用度值后，计算不同粒子的接收概率Pg:
[0082]
[0083] 其中，∑为求和运算符；
[0084] 在得到所有粒子的接收概率后，采用轮盘赌的方式获得群体最优位置xgbest。
[0085] 作为本发明所述的一种基于反射系数精确式的页岩VTI储层叠前混合反演方法进一步优化方案，步骤十一中，根据个体最优位置和群体最优位置来计算更新速度，进而采用迭代格式更新粒子包括：
[0086] 基于已获得个体最优位置 和群体最优位置xgbest，计算更新速度：
[0087]
[0088] 其中，vj,l+1、vj,l表示第j个粒子在第l+1次、第l次迭代中的更新速度，rand1(·)和rand2(·)表示两个介于[0,1]的随机数；κ2和κ3为两个控制更新步长的常数；xj,l表示第l次迭代的第j个粒子； 表示第l次迭代的第j个粒子的个体最优位置；
[0089] 基于更新速度，更新不同的粒子：
[0090] xj,l+1＝xj,l+κ1·vj,l+1
[0091] 其中，κ1为用以控制更新的步长的常数；xj,l+1表示第l+1次迭代的第j个粒子。
[0092] 本发明采用以上技术方案与现有技术相比，具有以下技术效果：
[0093] (1)针对目前页岩等VTI(具有垂直同相轴的横向各向同性)介质的反演现存的难点问题，提出了一种将线性和非线性优化算法相结合的混合反演策略，并将VTI介质的反射系数精确式引入作为正演算子，用以提高页岩等VTI介质的各向异性参数的反演精度；
[0094] (2)该混合反演策略将线性和非线性反演相结合，利用线性反演结果作为约束，生成较小的搜索范围和初始粒子，开展非线性反演以同时获取五参数的反演结果；该方法克服常规线性反演易落入局部极值的缺点，提高非线性反演的计算效率和准确性；采用反射系数精确式可克服近似表达式小角度、弱反射和弱各向异性的限制，提高反演的适用性。
附图说明
[0095] 图1是本发明混合反演算法的流程图。
[0096] 图2为基于VTI反射系数近似式的常规线性直接反演所获得的五参数结果。
[0097] 图3为本发明混合反演算法的线性反演结果。
[0098] 图4为基于线性反演结果所获得的搜索范围。
[0099] 图5为本发明混合反演算法的非线性最终结果。
具体实施方式
[0100] 下面结合附图对本发明的技术方案做进一步的详细说明：
[0101] 为克服常规线性反演易落入局部极值，而非线性反演搜索范围大计算效率低的问题，本发明提出了线性及非线性混合反演策略。为克服基于近似式反演远角度精度低、弱各向异性限制等问题，提出选用反射系数精确式的作为正演算子，以提高各向异性参数的反演精度。
[0102] 这里以一个测井模型测试为例进行说明：
[0103] 图1是本次申请一种基于VTI介质反射系数精确式的混合反演策略的实施流程，具体包括以下步骤：
[0104] 一、基于VTI介质反射系数精确式的正演
[0105] 在井口处获取弹性参数曲线，包括纵波速度α、横波速度β、密度ρ、各向异性参数ε和δ。同时设定该五参数为反演的目标参数，由五参数所组成的目标参数向量m为：
[0106] m＝[α1,...,αM,...1,...,βM,ρ1,...,ρM,ε1,...,εM,δ1,...,δM] (1)[0107] 其中，βi、βi、ρi、εi和δi(i＝1,....,M)表示为第i个时间采样点所对应的纵波速度、横波速度、密度、两个各向异性参数。
[0108] 从五个目标参数曲线中，通过Backus平均算法或者低通滤波来获取提取低频数据信息；而针对二维以及三维情况，采用井位置处的低频数据，并通过克里金插值法获得反演所用的三维初始模型。其中，地质层位解释被作为插值的横向约束。
[0109] 从PP地震角度道集的不同角度范围中提取不同角度下的地震子波序列，并利用子波序列建立子波矩阵W(θ)，表达式为：
[0110]
[0111] 其中，w(θi)表示θi角度所对应的子波矩阵：
[0112]
[0113] 基于褶积理论，合成地震记录dsyn可以利用子波矩阵W和反射系数序列r计算求得：
[0114] dsyn＝G(m)＝W·r (4)
[0115] G(m)为精确式正演算子的非线性实现过程，等同于合成地震记录。
[0116] 利用VTI介质反射系数精确方程可以计算反射系数r：
[0117] A·r＝b (5)
[0118] 其中，A和b均为中间过程矩阵，具有如下表达形式：
[0119]
[0120]
[0121] 其中，具有上标U和L的参数分别代表反射界面上层和下层的参数，而具有下标P和S的参数分别代表纵波和横波相关的参数。c11、c13、c33和c55为弹性刚度参数。符号p为水平慢度，s为垂直慢度。 n均为中间变量参数，具体计算公式如下：
[0122]
[0123]
[0124] 其中，
[0125]
[0126] 式中，sP和sS表示纵波和横波所对应的水平慢度。
[0127] 上面反射系数精确式均为刚度参数c11、c13、c33、c55的函数，而纵波速度α、横波速度β、密度ρ，各向异性参数ε和δ曲线，与刚度参数c11、c13、c33、c55之间的关系如下：
[0128] c33＝α2ρ,c55＝β2ρ (11)
[0129] c11＝(2ε+1)α2ρ,c13＝ρη‑β2ρ (12)
[0130] 其中，η为中间变量
[0131]
[0132] 二、基于VTI反射系数精确式的叠前线性反演方法
[0133] 建立反演目标函数S(m)，表示为：
[0134]
[0135] 其中，d为实际叠前地震道集数据，G(m)为精确式正演算子的非线性实现过程，即合成地震记录。λ为正则化参数，Cm为目标参数向量计算得到的协方差矩阵，u为目标参数向量的期望。
[0136] 计算反射系数精确式正演算子G(m)对参数向量m的偏导数
[0137]
[0138] 其中，W为子波矩阵， 为反射系数对目标参数向量的偏导数，包括对反射界面上U L层介质模型参数m和下层介质模型参数m的偏导数 和 具体的计算公式为：
[0139]
[0140]
[0141] 式中，A和b为求取反射系数r时的中间过程矩阵。 分别U L
表示中间变量矩阵A和b对反射界面上层介质模型参数m和下层介质模型参数m的偏导数。
[0142] 和 的具体计算表达式为：
[0143]
[0144] 式中，上标U和L表示反射界面上层介质和下层介质所对应的参数。U
和 分别为中间变量 和 对上层介质模型参数m的偏导数。
L U U U U
和 分别为中间变量 和 对下层介质模型参数m 的偏导数。a 、b 、d 、e的
具体计算式如下：
[0145]
[0146]
[0147]
[0148]U
[0149] 式中， 分别表示上层介质慢度参数 p对模型参数m的U L
偏导数。 为上层介质刚度参数 对模型参数m的偏导数。a 、
L L L U U U U
b、d、e与a、b、d、e具有同样的计算形式。
[0150] 和 具有一样的计算形式，这里仅展示 计算表达式：
[0151]
[0152] 式中， 分别为中间变量 对上层介质模型参数mU的偏导数。fU的具体表达式如下：
[0153]
[0154] 计算目标函数S(m)的雅可比矩阵J(m)，即一阶偏导数矩阵：
[0155]
[0156] 其中， 为精确式正演算子对目标参数向量m的偏导数。
[0157] 针对目标函数S(m)的伪海森矩阵H(m)，采用了内部小循环迭代形式来计算获得：
[0158] 对于线性反演循环迭代中第k次迭代、内部小循环首次迭代的伪海森矩阵H k,1(m)，可由下式计算得到：
[0159]
[0160] 对于线性反演循环迭代中第k次迭代、内部小循环第i次迭代的为海森矩阵Hk,i(m)，可由下式计算获得：
[0161] Hk,i(m)‑1＝(Qk,i‑1)THk,i‑1(m)‑1(Qk,i‑1)+Kk,i‑1sk,i‑1(sk,i‑1)T (27)[0162] 其中，Hk,i‑1(m)为外部大循环第k次迭代内部小循环第i‑1次迭代的伪海森矩阵；外部大循环是指线性反演循环，Qk,i‑1、Kk,i‑1、sk,i‑1为外部大循环第k次迭代内部小循环第i‑1次迭代的中间过程矩阵，计算形式如下：
[0163] Qk,i‑1＝I‑Kk,i‑1yk,i‑1(sk,i‑1)T (28)[0164] Kk,i‑1＝1/(yk,i‑1)T sk,i‑1 (29)
[0165] 式(27)至(29)中，
[0166] sk,i‑1＝mk+1‑mk，yk,i‑1＝J(mk+1)‑J(mk) (30)
[0167] 上式中，mk为线性反演循环迭代中第k次迭代的目标参数向量，J(mk)为线性反演循环迭代中第k次迭代的雅可比矩阵。
[0168] 线性反演循环迭代的第k次迭代时，目标参数向量的更新方向‑gk(m)由下式计算得到：
[0169] ‑gk(m)＝‑Hk(m)·Jk(m) (31)
[0170] 其中，Hk(m)和Jk(m)分别代表外部大循环中第k次迭代所对应的伪海森矩阵和雅可比矩阵。
[0171] 利用如下的更新迭代格式来更新参数向量mk+1：
[0172]
[0173] 其中，mk+1、mk分别表示外部大循环中第k+1次、第k次迭代的目标参数向量，vk、vk+1分别表示外部大循环中第k次、第k+1次迭代的中间过程矩阵，是动量系数，a为更新步长。
[0174] 若线性反演循环迭代的迭代次数未达到最大迭代次数，则开始下一次反演迭代，重复上述内容；若达到，则结束线性反演迭代过程，输出线性反演结果
[0175] 三、基于VTI反射系数精确式的叠前非线性反演方法
[0176] 基于线性反演结果以生成较小范围的搜索区间。以模型中某一深度样点为例，用Li min maxm 表示线性反演结果，则基于线性反演结果可生成搜索范围[m m ]:
[0177] 当|mRe‑mLi|＜b·mLi (33)Re Li Li
[0178] 当|m ‑m |≥b·m (34)
[0179] 其中，mmin为搜索范围的最小值，mmax为搜索范围的最大值，mRe表示井口处抽取到的真实参考值，b为尺度因子，可随参数不同而变化。
[0180] 基于线性的反演结果设定初始粒子。设定初始粒子时，粒子群x可以表示为：
[0181] x＝[x1,...,xj,...,xN] (35)
[0182] 其中，xj表示第j个粒子，也是第j个可能的模型参数；并将线性反演结果mLi赋值给N个粒子中的任意一个，最终形成初始粒子群。
[0183] 将第j个粒子xj代入目标函数计算得到S(xj)，S(xj)也被称为是第j个粒子的适用度；计算出所有粒子的适用度后，选择适用度最小的粒子赋值为个体最优位置
[0184] 当k＝1,
[0185] 当k>1,
[0186] 其中，l为非线性反演循环迭代次数，xj,l表示第l次迭代中的第j个粒子，S(xj,l)、S(xj,l‑1)分别为第l次、第l‑1次迭代中的第j个粒子的适应度。
[0187] 计算替换适用度Zr(xj,xgbest,Tl)，用以计算接收概率Pg：
[0188]
[0189]
[0190] 其中，xj表示第j个粒子，xgbest为群体最优位置，S(xj)和S(xgbest)分别表示将xj和gbest gbestx 代入目标函数计算得到的适用度； 为xj和x 适用度的差，l为非线性反演循环迭代的迭代次数，Tl为第l次迭代中的温度，τ表示为一个常数；
[0191] 计算好所有粒子的替换适用度后，可计算获得不同粒子的接收概率Pg:
[0192]
[0193] 其中，∑为求和运算符。在得到所有粒子的接收概率后，采用轮盘赌的方式获得群gbest gbest体最优位置x 。基于已获得个体最优位置 和群体最优位置x ，计算更新速度：
[0194]
[0195] 其中，vj,l+1，vj,l表示第j个粒子在第l次迭代中的更新速度，rand1(·)和rand2(·)表示两个介于[0,1]的随机数；κ2和κ3为两个控制更新步长的常数；xj,l表示第l次迭代的第j个粒子； 表示第l次迭代的第j个粒子的个体最优位置。
[0196] 基于更新速度，更新不同的粒子：
[0197] xj,l+1＝xj,l+κ1·vj,l+1 (42)
[0198] 其中，κ1为用以控制更新的步长的常数；xj,l+1表示第l+1次迭代的第j个粒子。若迭代次数未达到最大迭代次数，则开始下一次反演迭代，重复上述内容；若达到最大迭代次数，则结束非线性反演，输出所有粒子，选取其中最好结果作为最终的非线性反演结果。最好结果就是粒子所对应的适应度值最低的结果。
[0199] 图2为基于VTI反射系数近似式的线性直接反演所获得五参数结果，包括纵波速度α、横波速度β、密度ρ、各向异性参数ε和各向异性参数δ；其中，黑色实线为真实测井模型，黑色虚线为初始模型，灰色点线为反演结果；可以看出近似式直接反演获得结果并不理想。图3是本发明混合反演算法线性反演部分获得结果，同样包括有α、β、ρ、ε和δ；其中，黑色实线为真实测井模型，黑色虚线为初始模型，灰色点线为反演结果；可以看出该结果比图2中近似式结果有很大提升，但是与实际模型相比仍有很大误差。图4为利用线性反演结果生成的搜索范围；其中，黑色实线为真实测井模型，灰色点线为线性反演结果，灰色虚线为生成的搜索范围；可见由于线性结果大程度贴近真实值，基于此可以生成较小范围的搜索区域，从而提升后续非线性反演的效率与精度。图5为本发明混合反演方法的非线性反演最终结果；
同样地，黑色实线为真实模型，灰色点线为反演结果；相比于图3中线性反演结果，非线性结果有很大的提升，与真实模型之间的一致性更好。
[0200] 以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围内。
